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the Estimation of Laminar Skin Friction, 
including the Effects of Distributed Suction 


N. CURLE, M.Sc., Ph.D. 


(Aerodynamics Division, National Physical Laboratory) 


SUMMARY: Stratford's analysis of the laminar boundary layer near separation 


uses two physical ideas. In the outer part of the boundary layer, where viscous 
effects are small, the development is given by the condition that the total head 
is constant along streamlines, apart from a second-order correction for 
viscosity. Near the wall, however, viscous forces must balance the pressure 
forces, and the profile adjusts itself accordingly. Quantitatively these ideas 
yield a simple formula for predicting separation, which has been found to be 
particularly accurate 


In this paper it is indicated how the same approach may be used to yield 
the full distribution of skin friction along the wall. Further, the etfects of 
suction may be incorporated into the method. Physically, suction atiects the 
outer part of the boundary layer in that the streamlines are drawn towards 
the wall when suction is applied. At the wall, the balance between viscous and 
pressure forces 1s influenced by the momentum of the fluid which is 
awat When these eflects are accounted for quantitatively, the 
formula for the skin friction is sull very simple 


sucked 


rp il 
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Several examples are considered, and comparison is made with exact 


theory and with approximate results by other methods. It is indicated that 
the method has a useful range of validity 


1. Introduction 


Flow in the laminar boundary layer on a solid body has been treated by 
numerous workers, and many methods of solution have been devised, of varying 
speed and accuracy. (For a detailed summary, see Ref. 1.) If only the position of 
separation is required, the best method for obtaining accurate results is probably 
that due to Stratford If, in addition, it is necessary to estimate the detailed 
distribution of skin friction or the momentum and displacement thicknesses, then 
the method of Thwaites™’ is very eflicient. The numerical accuracy of Thwaites’s 
method and the speed of Stratford’s method were both later improved by Curle and 
Skan'”. If, however, one wishes to know accurately the shape of the boundary- 
layer profile, as is necessary for a stability calculation, then a more powerful (and 
more lengthy) method must be used, such as that due to Head’ 


When one considers the boundary layer with suction (or blowing), accurate 


methods are hard to find and the only method in which one can really place any 
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The Estimation of Laminar Skin Friction, 
including the Effects of Distributed Suction 


N. CURLE, M.Sc., Ph.D. 


(Aerodynamics Division, National Physical Laboratory) 


Summary: Stratford's analysis of the laminar boundary layer near separation 
uses two physical ideas. In the outer part of the boundary layer, where viscous 
effects are small, the development is given by the condition that the total head 
is constant along streamlines, apart from a second-order correction for 
viscosity. Near the wall, however, viscous forces must balance the pressure 
forces, and the profile adjusts itself accordingly. Quantitatively these ideas 
yield a simple formula for predicting separation, which has been found to be 
particularly accurate. 


In this paper it is indicated how the same approach may be used to yield 
the full distribution of skin friction along the wall. Further, the effects of 
suction may be incorporated into the method. Physically, suction affects the 
outer part of the boundary layer in that the streamlines are drawn towards 
the wall when suction is applied. At the wall, the balance between viscous and 
pressure forces is influenced by the momentum of the fluid which is sucked 
away. When these effects are accounted for quantitatively, the resulting 
formula for the skin friction is still very simple. 


Several examples are considered, and comparison is made with exact 
theory and with approximate results by other methods. It is indicated that 
the method has a useful range of validity. 


1. Introduction 


Flow in the laminar boundary layer on a solid body has been treated by 
numerous workers, and many methods of solution have been devised, of varying 
speed and accuracy. (For a detailed summary, see Ref. 1.) If only the position of 
separation is required, the best method for obtaining accurate results is probably 
that due to Stratford. If, in addition, it is necessary to estimate the detailed 
distribution of skin friction or the momentum and displacement thicknesses, then 
the method of Thwaites® is very efficient. The numerical accuracy of Thwaites’s 
method and the speed of Stratford’s method were both later improved by Curle and 
Skan™. If, however, one wishes to know accurately the shape of the boundary- 
layer profile, as is necessary for a stability calculation, then a more powerful (and 
more lengthy) method must be used, such as that due to Head’. 


When one considers the boundary layer with suction (or blowing), accurate 
methods are hard to find and the only method in which one can really place any 
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confidence is that of Head’. As Head points out, any one-parameter method, 
which may give quite reasonable answers for solid boundaries, is at a distinct 
disadvantage when suction is applied, whereas his two-parameter method remains 
almost as accurate. The reason for this is that, if finite suction is applied at some 
position, then the curvature of the velocity profile at the wall will change accord- 
ingly, although the momentum thickness, for example, will be continuous. These 
independent conditions cannot, in general, be satisfied by a one-parameter family 
of profiles, for which reason the methods of Thwaites® and Schlichting’s'” 
adaptation of the classical Pohlhausen™ method give only a rough indication of the 
effects of suction. 


The purpose of this present paper is twofold. Firstly it is indicated how the 
method used by Stratford can be extended so as to yield not only accurate 
predictions of separation, but also the detailed distribution of skin friction with 
equal accuracy. Secondly, and simultaneously, it is shown how the effects of 
distributed suction may be incorporated into the theory. Stratford’s method was 
originally designed for dealing with flows in which a large pressure gradient caused 
a sharp curvature of the profile near the wall, while not greatly affecting the general 
shape in the outer part of the layer. One might reasonably hope, therefore, that 
the method would deal with sharp local profile curvatures caused by suction and 
so lead to a particularly speedy method. 


The generalisation of Stratford’s theory for retarded flow is first developed 


(Section 2), including the allowance to be made for an initial favourable pressure 
gradient (Section 3). Three distinct types of example are then worked out, the 
solutions in each case being compared with exact solutions or with accurate solutions 
by other approximate methods. 


For flows without suction (Section 4) it is known that Stratford’s method will 
usually yield accurate predictions of separation.“’ The detailed distribution of 
skin friction is derived in this paper for the problem in which the main stream 
velocity u, decreases linearly with distance x along the boundary“, so that 


u, =u, (1—-x/D, : ‘ (1) 


where u, and / are a representative velocity and length respectively. In addition 
the skin friction is predicted in the region of adverse pressure gradient for 
Schubauer’s experimental velocity distribution’. In each case the agreement with 
the exact solution is good. Further, a case with a uniform favourable pressure 
gradient is considered, corresponding to a main stream velocity 


No exact solution is available for this problem, but, as Thwaites’s method usually 
yields accurate predictions of skin friction in regions of accelerated flow, a solution 
by his method has accordingly been obtained for comparison, and the two solutions 
give values for the local skin friction which differ nowhere by more than 7 per cent. 
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LAMINAR SKIN FRICTION WITH SUCTION 


Further, the solution for large x may be compared with the exact asymptotic 
solution, corresponding to one of the Falkner-Skan similarity profiles’*’. The 
predicted skin friction is in error by only 3 per cent. 


The method is then used (Section 5) for predicting separation in problems with 
prescribed suction and pressure gradient. We consider the linearly retarded main 
stream (equation (1)), with a uniform velocity of suction w,. The position of 
Separation is estimated for various values of a suction parameter o =w, [I//(u,¥)]'/?, 
where v is the kinematic viscosity. The results are compared with two solutions 
by approximate methods due to Thwaites: ”, and appear to be quite reasonable. 
Further, an estimate is made of the suction distribution which would just prevent 
separation, the main stream velocity (equation (1)) again being used. No exact 
solution to this problem is available, but a solution by Head“, which is almost 
certainly quite accurate, is compared with the prediction of this paper. The present 
method indicates greater suction requirements immediately downstream of where 
suction is first applied, but the agreement is better some distance downstream. 


The method is then applied (Section 6) to flows with zero pressure gradient 
and prescribed suction, the distribution of skin friction being predicted. The case 
of uniform suction is considered, in which the profile, of the Blasius type at the 
leading edge, gradually deforms into the asymptotic suction profile far downstream. 
The solution is compared with the exact solution of Iglisch’”. For values of 
w, [2x/(u,»)]'/? less than about 0-2, the predicted skin friction is in error by less 
than 5 per cent, at a position where it is about 30 per cent greater than in the absence 
of suction. The more general problem is then considered, with a solid entry length 
and uniform suction applied downstream of a prescribed position. Here the Blasius 
solution is valid up to the position where suction is applied, and deforms downstream 
of this position into the asymptotic suction profile. 


For flow without suction the method is thus both accurate and speedy. For 
flows with suction, however, the shortage of accurate numerical solutions makes it 
difficult to assess accuracy. What evidence is available indicates that the method 
should in general be at least as accurate as most of the one-parameter methods, as 
well as being considerably quicker to use. The question of accuracy is discussed at 
various points in the text. 


NOTATION 

x,z Cartesian co-ordinates, parallel and perpendicular to wall 

u,w velocity components parallel to x- and z-axes respectively 

1,1, representative lengths 

o,é suction parameters, defined in equations (46) and (49) respectively 
p pressure 
viscosity 
p density 
v kinematic viscosity, 
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stream function 
w, Velocity of suction, = —w (0) 

7 skin friction at wall 

a_ a function introduced in equation (10) 

n,k constants 

C, pressure coefficient, =(p — p,)/(4pu,”) 
xX, a length, defined by equation (30) 


Suffixes 
0 values at leading edge or at start of pressure rise 
1 local values in main stream 
w_ values at the wall 
B_ Blasius profile 
j position of join of inner and outer profile 


m_ values at pressure minimum 


2. Development of General Theory for Retarded Flow 


The inner and outer portions of the boundary layer are considered separately. 
In deriving the theory attention will first be restricted, following Stratford, to flows 
in which the pressure is constant between x=0 and x= x,, with a sharp but otherwise 
arbitrary pressure rise starting at x=x,. In addition, we shall allow an arbitrary 
distribution of suction when x > x,. If the pressure gradient is very great, then 
viscous forces will be dominated by pressure forces, except in a very thin layer near 
the wall. Hence it is permissible to consider the development of the outer layer by 
first assuming inviscid flow and then adding on a secondary correction for the effects 
of viscosity. 


If we assume nearly inviscid flow when x > x,, z =z), where z is the distance 
measured in a direction normal to the wall and suffix “ j” refers to values at the join 
of the inner and outer layers, then in this region the total head is nearly constant 
along streamlines. 


Thus —(p — py) + AH, ‘ (3) 
where the stream function vy is given by 


z 
Y= w, dx+ | u dz, . (4) 
Xo 4 
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u is the velocity parallel to the wall, p is the density and p the pressure of the fluid, 
w, is the velocity of suction and AH is the change in total head due to viscosity. 
Now the shape of the velocity profile in the outer layer is not greatly changed by a 
sharp pressure rise of the type being considered or by suction. Consequently the 
viscous forces are given approximately as if there were no pressure rise or suction. 
Accordingly the change of total head due to viscosity is 


AH = (4puy’)x.. (5) 


where u, is the Blasius solution for a boundary layer with no pressure gradient. If 
this is superimposed on equation (3), noting that u=u, when x=.x,, then the outer 
layer is given by 

(4pu’)x.. =(4pug*)x.y —(p Py). 


which may be written 


u? (x, = up? (x, — 


This result may also be deduced algebraically from the boundary-layer equations, 
provided that terms of order (x—x,)’ may be neglected, as Stratford showed. 
Further, by similar arguments, which amount essentially to the differentiation of 
equation (6) with respect to ¥, it may be deduced that in the outer layer 


(=) 


u 
027/ Uy O27 /, 


In the inner layer, in contrast, viscous forces become very important, so that 
when there is no suction the resultant viscous force on an element of fluid at the wall 
must exactly balance the pressure-gradient force. When there is suction, momentum 
is being transferred across the boundary and the equation representing the balance 
of forces at the wall is modified. This equation, which may be obtained by letting 
z—> 0 in the momentum equation, is 


ou oy dp 
f AZ) ~ 


It is assumed, generalising Stratford’s analysis, that the velocity profile in the inner 
layer may be adequately represented by a form 


uu We) 28+ az" = 0) 
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Here - is the skin friction » (@u/0z), and, for the moment, a and the constant n are 
arbitrary. The condition (9) is, of course, identically satisfied by (10), which holds 
for x > X,, ZS 


It is now postulated that ¥, u, Cu/dz and 0?u/0z* shall be continuous at the 
join of the inner and outer layers. These four conditions are sufficient to determine 
7, a, Zz, and v,, with n appearing as an arbitrary parameter. As in Stratford’s theory, 
the simplifying assumption is made that the inner layer is sufficiently thin for the 
Blasius profile to be linear throughout it. This will be true provided that 


Uy (x, 
u 


where k is, say, 0-6, and must certainly be less than unity. The implications of this 
restriction will be examined later. Assuming that it holds, the Blasius profile is 
approximated by 


Ir 1/2 
_ TR 
od (57) 2. 


where 7, is the skin friction for the Blasius layer, and relations (12) to (14) hold for 
Y<y,. The joining conditions then become 


= — Po) (15) 
iL 
uv ufos dx+4;z dp 24 (16) 
1 TW, 
dx 
dx 


The solution of these five equations is fairly straightforward, and leads to the 
following result. 
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IC 
Votes 2 
fc 2: 
x { 4 {c,- Ww, dx} 
dx uu,” 


(n—1) (n?+3) (1- {1 3 (n? —1) 


0:00405 (n — 1)° (n? + 3) 3 (n?—1) 4 


n* (n+ 1)(n—2)° 


where C, is the pressure coefficient, defined as (p — p,)/(4pu,”). 


A criterion for separation is obtained by setting = =0 in equation (20). If this 
is done, then for the case of solid boundaries, w,=0, it follows that 


0:00405 (n — 1)* (n? + 3) 
Lx n? (n+ 1)(n—2) (21) 
When n=6 the right-hand side becomes 0:0049, and when n= 4 it becomes 0-0065, 
the values obtained by Stratford. Further, when n=3 it becomes 00-0108. Now 
Curle and Skan™ showed that the best all-round accuracy in predicting separation 
would be obtained by taking the right-hand side of equation (21) equal to 0-0104, 
which corresponds to 


n=3-043. . ‘ - 


Accordingly the theory will be developed assuming this value of nm, so that the 
equation (20) for - becomes 


_ 


2 
dx 


uu,” J 


The method is nominally applicable, according to condition (11), when 


Ug (x, ¥;) 
Uy 


and, by relation (12), this condition becomes 


February 1960 7 


Xp 


N. CURLE 


where it can be shown that 


x 


{ w, dx. 


, 916(1+1-017) ax} + 2% 


When there is no suction, the method should be valid at separation provided that 


Since k is about 0:6, and certainly is less than unity, it follows that the method 
should strictly apply only when (C,)., is less than about 0:04, with an absolute 
upper limit of about 0-11, and for greater values of C, it should be necessary to use 
the exact joining conditions rather than the simpler ones which assume the linearity 
of the Blasius profile. Nevertheless, it has been found that the simpler method 
yields accurate predictions of separation for values of C, as great as 0:25. It will 
be found in this paper that the method also applies far beyond its nominal range in 
predictions of the skin friction 7, and for flows with suction. 


3. Allowance for Initial Favourable Pressure Gradient 


The theory as developed in Section 2, although it does not exclude all 
accelerated flows, does assume that the initial velocity profile is of the Blasius type, 
which will be true if the leading edge is sharp. Boundary layers beginning at a 
stagnation point require different treatment. 


Consider, then, a flow in which the main stream velocity u, (x) increases from 
zero at a stagnation point x=0, reaching a maximum (at the pressure minimum) and 
then decreasing again. Stratford indicated that at the position of the pressure 
minimum, x=x,,, the velocity profile would be very similar to the Blasius profile of 
an equivalent distance x,,—x,. The position x= x, of the equivalent origin is deter- 
mined by equating the known momentum thickness of the equivalent Blasius profile 
to the true boundary-layer momentum thickness at x= X,,. 


Now, according to Thwaites’s method, which predicts momentum thickness 
accurately, the true momentum thickness at a position x is given by 


x 


6,2=0 45vu,~* | dx, . 


from which it may be shown that 


x 


ax. 


m 
Um 
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The development of the boundary layer in x > x,, will be given by the following 
modified form of equation (23) :— 


x 


(x - = = T \ ( w, dx =0-0104 (1 —T) (1 + 2-027), 
dx } J 


Xm 


(31) 


Tp 0:33206 


_= (32 


where 


and C, is also defined in terms of the pressure and velocity at x= x,,. This form of 
the equation will be used in the particular application to flow past an elliptic 
cylinder (Section 4.2). It is to be noted that, for flows where the boundary layer 
begins at a stagnation point, the method cannot deal with suction upstream of the 
pressure minimum. This is not a serious drawback, from a practical point of view, 
as suction would normally be applied only in a region of retarded flow. 


4. Application to Some Flows Without Suction 


When there is no suction the general formula (31) becomes 


(x— x1)" .  . (3) 


dx 


In this section three different velocity distributions are considered, and in each case 
T (x) is determined for comparison with exact theory. 


4.1. LINEARLY RETARDED VELOCITY 


In this example the main stream velocity is 
u,=u,(1—x/) 


Since the flow is everywhere retarded, it follows that x, = 0. 
Also C, 


so equation (33) becomes 


(1- (2- 0-0026 (1 —T)' (1+ 2-027). 


The values of T have been computed for several values of x//, and values of 
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Ficure 1. Skin friction for linearly retarded velocity distribution. 


have been deduced therefrom. The results are shown in Table I and Fig. 1, as are the 
exact values obtained by Howarth". The agreement, it will be noted, is very good, 
the distribution of skin friction being given accurately, except just upstream of the 
separation point, which, none the less, is determined to an accuracy of one per cent. 


The values given by the approximate method of Thwaites’, generally agreed 
to be one of the more accurate approximate methods, are of roughly the same 
accuracy as the present results. 


4.2. SCHUBAUER’S EXPERIMENTAL DISTRIBUTION 


Hartree“? has examined the experimental pressure distribution obtained by 
Schubauer for flow past an elliptic cylinder. He has carried out a precise numerical 
integration of the boundary-layer equations, assuming Schubauer’s experimentally 
obtained values of C,, and prescribing dC,/dx by numerical differentiation. With 
the minor axis of the ellipse as the representative length /, the pressure minimum is 
at x/1=1-300, and the velocity maximum is u,,=1-295u.,, where u,, is the velocity 
at infinity. The position of separation is x//=1-983. The change of origin x,, to 
account for the initial favourable pressure gradient, is easily calculated to be 
x; =0-3941, whereupon the value of T can be obtained from equation (33) at each 
station x// => 1-3 for which C, is given. From this is calculated 


2 (1-295) G) 9791 (* 


which is equal to Hartree’s “20*y/¢y*” in his Table 8. The results are shown in 
Table II, and again the agreement is extremely good. The skin friction is predicted 
to within 1 per cent for 1-3 <= x/1 <= 1-8, and the position of separation is also given 
to an accuracy of 1 per cent. 


The accuracy in this case is rather better than that obtained by Thwaites’s 
method. In general one would expect this to be so, since Stratford’s method 
predicts the separation position even more accurately than Thwaites’s method". 
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0 
00125 
0:0250 
0:0375 
0:0625 
0:0750 
0-0875 
0°1000 
0°1125 
0-1200 


TABLE | 


SKIN FRICTION FOR LINEARLY RETARDED VELOCITY 


(SECTION 4.1) 


( | €) 
Cz 


Predicted Value Exact Value Thwaites® 
x x 
2°733 2°739 
1-762 1-772 1-765 
1-295 1-309 
0-995 1-011 1-019 
0-772 0°790 
0-613 0631 
0-435 0-459 
0-289 0-315 
0-140 0-163 
0-029 0-000 0-072 
TABLE II 


SKIN FRICTION FOR SCHUBAUER’S PRESSURE 


SKIN FRICTION FOR UNIFORM FAVOURABLE P 


Values of T 


(SECTION 4,2) 


2 (1-295)3,2 ( ) (=) 
Un Cz w 


Predicted Value Exact Value 


0-958 0-943 
0°829 

0-685 0-687 

0°399 0-395 
0:314 

0:230 

0-143 

0-033 


TABLE III 


(SECTION 4,3) 


DISTRIBUTION 


Thwaites‘®) 


0-912 


0°346 


RESSURE GRADIENT 


Values of T 


Equation (39) 


2°718 
3-124 
3-318 
3-509 
3-697 


| 
Equation (39) Thwaite® | x/l 
| 
1-000 1-000 0:7 
1-298 1-263 | 0-8 
1-568 1-519 0-9 
1-820 1735 «| 10 
2-058 91s | 
2-287 2-133 | 1-2 
2-507 2-377 | 
11 


2°810 
3-004 
3194 
3-378 
3-561 


1:30 
1-40 
1-50 
1°60 
1:70 
1:80 
1:85 
190 
1:95 
2:00 
— 
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4.3. UNIFORM FAVOURABLE PRESSURE GRADIENT 


We consider now the accelerated velocity distribution 


which corresponds to a pressure coefficient 


Although the flow is everywhere accelerated, x, =0, nevertheless, since there is no 
stagnation point. So equation (33) becomes 


(x/1)’ =0:0104 (T —1)° (1+2-027), . (39) 


from which T may be deduced as a function of x//, as shown in Table III. No 
exact solution is available for purposes of comparison, but it may be deduced from 
equation (39) that, for small x//, T may be expanded in a series 


and that, for large x//, the asymptotic form is 


Fuld ‘ . (41) 


Now the exact values of the numerical coefficients in relations (40) and (41) may be 
found as follows. The coefficient of x// in (40) is found from Howarth’s results, 
and is 3:07. The coefficient of (x//)°'* in (41) is deduced from the Falkner-Skan 
similarity solutions, and is 2:705. It would therefore seem a reasonable conclusion 
that the skin friction predicted by the present method is nowhere in error by more 
than about 3 per cent. This conclusion is reinforced by comparison with a solution 
which the author has obtained by Thwaites’s method“, usually accurate for 
accelerated flow. If the Thwaites solution is expanded as in (40), then the 
coefficient of x// is 2:91, so that for small x the error in the present solution should 
be about five-eighths of the error in Thwaites’s solution. As shown in Table III, 
the greatest discrepancy between the two solutions is about 7:5 per cent; one would 
conclude that at worst ine predicted skin friction becomes 3 per cent too high before 
tending to its asymptotic valve of 3 per cent too low. 


It is worthy of note that, in each of the three examples considered in this 
section, the method applies far beyond its nominal range of application. In the 
first two examples the pressure coefficient at separation is about 0°23 and 0-09 
respectively, compared with the nominal limit of about 0:04. In the third example, 
the function which should be less than k*? ~ 0-36 increases from about 1-87 at the 
leading edge to about 1-70 (x//)''* for large x//. One can only conclude that in 
some fortuitous manner the errors caused by the various approximations made in 
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the theory, namely the approximation to the outer profile and the assumption of a 
linear Blasius régime throughout the inner profile, must cancel out, at any rate as 
regards the prediction of skin friction. 


5. Predictions of Separation for Some Flows with Suction 
5.1. THE SIMPLIFIED EQUATION 
This section is concerned only with separation properties, the assumption being 


made that separation occurs where the skin friction is zero. The general equation 
(23) accordingly simplifies to 


x 


(&) {c, j Ww, dx} 0:0104, (42) 


provided that the velocity distribution is everywhere retarded. Further, we consider 
only cases of a linearly retarded velocity 

u, =u, (1—x/D), ‘ ‘ ; . (43) 
for which C, is given in equation (35). Hence equation (42) becomes 


x\? x x 066412 ( } 
(1 w, dx ):0026 (44) 


5.2. UNIFORM SUCTION 


When uniform suction w, is applied from the leading edge, then equation (44) 
indicates that separation occurs where 


"(1 *) {27 -0-664127(*) } -0:0026, . (45) 


where « is the suction parameter 


1/2 
) : . (46) 


The solution of equation (46), giving the relationship between the suction 
parameter « and the position of separation, is shown in Table IV and Fig. 2. It may 
be noted from Table IV that as o increases the position of separation moves back, 
from x//=0:12 when =0 to x/1=0-18 when =0°74, and then to x//=0-22 when 
=0:974 and to x/1=0-30 when «=1-240. Beyond this point the position of 
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FiGuRE 2. Position of separation for linearly retarded velocity, with uniform suction. 


separation moves back very rapidly as ~ increases, so that an increase in o from 
1:240 to 1°545 results in separation being delayed from 0:30 to 0-60. When 
o=1:5466 separation just occurs at x//=0-62, and separation is completely 
suppressed when o > 1:5466. 


The result that separation is completely suppressed for large enough @ is not 
new. For example, Thwaites’> suggested a critical value o=4, and earlier 
Prandtl"”? and Preston” had suggested 2:18 and 1°61 respectively. The correct 
value of « for complete suppression of separation will, of course, depend upon the 
pressure distribution. The possibility that Thwaites’s value may be too high can 
be reasonably inferred from some of the results of Head. In Fig. 10 of Ref. 13 
Head shows the results of his method and Thwaites’s method for the problem (to be 
considered in Section 5.3) of what suction distribution would be required to produce 
a zero-skin-friction layer with a linearly retarded velocity distribution. Thwaites’s 
curve of v,*, which is just the local value of «, is tending to the limiting value 4 as 
x/1—> 1, as would be expected. On the other hand, Head’s values are smaller 
when x// > 0-16, and appear to be tending to a limit of about 2. This suggests 
that possibly the critical value « ~ 1:55, as given by the present method, may not 
be too wildly in error. 


It is interesting to note the manner in which suppression of separation is 
predicted to arise. As o increases, separation moves back to x/]=0°62 when 
o =1-5466, and no separation occurs if « is further increased; separation in 
0-62 < x/1< 1-00 cannot be achieved with uniform suction. Whether or not this 
really happens will presumably only be decided when an exact solution is obtained. 


In Fig. 2 the predictions of the present theory are compared with the results of 
two theories due to Thwaites‘ ”. It is seen that the present results agree reasonably 
with both when « < 1-2, so presumably they are accurate in this range of o at least. 
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TABLE IV 


POSITION OF SEPARATION FOR LINEARLY RETARDED VELOCITY, WITH UNIFORM SUCTION 


(SECTION 5.2) 


o 


0 
0-326 
0°739 
0°870 
0:974 


TABLE V 
SUCTION REQUIRED FOR ZERO SKIN FRICTION LAYER 


(SECTION 5.3) 


x 


( wdx w, 


Equation (44) Head\'*) Equation (44) Head") 


— 0-004 0-000 
0°65 

0:046 0-013 
0°88 

0-091 0:030 


0°133 0-051 


0:174 0-076 
0:273 
0°372 


5.3. SUCTION REQUIRED TO PRESERVE ZERO SKIN FRICTION 

In the absence of suction the boundary layer with external velocity (equation 
(43)) separates where x//=0-120, the solution of equation (44) being x//=0-121 
when w,=0. If the correct amount of suction is applied downstream of this position 
it is possible to maintain zero skin friction. The required suction velocity w, is given 
as a function of x by equation (44). 


By solving this equation, the value of (u,v/)~*/ ( w. dx was obtained as a 

function of x/1, and then w, [//(u,v)]'/? by differentiation. The results are shown in 
Table V and Figs. 3 and 4, where they are compared with the results obtained 


by Head"®’. 
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= x/l | x/l 
0-121 1:240 0-30 
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0:18 1545 0-60 
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0:22 
x/l 
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0-18 
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Ficure 3. Suction quantity for zero skin friction with linearly retarded main stream. 


It is noted that the present method agrees with Head’s in predicting that a 
non-zero suction must be applied at the position where separation would otherwise 
occur, namely x//~0:12. The magnitude, however, is overestimated by the 
present method, by a factor of 4. Downstream, however, the results in Fig. 4 are 
in much better agreement. For example, at x// ~ 0-17, the predicted suction is 
twice Head’s value, and, by extrapolating Head’s results it is estimated that the 
predicted suction is only 25 per cent too high when x// =~ 0:25 and about 7 per cent 
too high when x// 0:30. 


The results for the suction quantity | w, dx behave appropriately, the curves 
0 
in Fig. 3 being roughly parallel when x// > 0-25. 


It is not absolutely certain just how accurate Head's results are, but it is quite 
certain that, in the absence of an extremely lengthy numerical integration of the full 
boundary-layer equations, they constitute the most reliable solution available. 
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Further, it seems almost certain that the required suction velocity w, should 
increase with distance, as the following argument suggests. If uniform suction were 
applied in x >0-12/, then presumably separation would still occur at some down- 
Stream position where greater suction would be required if separation is to be 
further delayed. This argument is not conclusive, of course, since the magnitude of 
the adverse pressure gradient is decreasing in the present example, but it does 
strengthen the case for accepting Head’s results rather than those of the present 
paper in the region 0-12 <= x/1/<0-25, where there is a discrepancy between them. 
Accordingly it is concluded that the present method does not yield accurate pre- 
dictions of suction near the position where suction begins, for this example at any 
rate, but does so farther downstream. 


A possible explanation of the initial discrepancy lies in the fact that the inner 
profile is rather thick in this example, so that the quantity of equation (27), which 
should be less than (0°6)* if the method is theoretically applicable at separation, is 
equal to 2:01 at the initial point of zero skin friction. Now Head points out that 
the main drawback of one-parameter methods lies in their implicit assumption that 
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Ficure 4. Suction velocity for zero skin friction with linearly retarded main stream. 
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changes in suction velocity affect the whole boundary layer, whereas in fact the 
effects can only be felt as they diffuse outwards. Similarly, in the present method 
it is implicitly assumed that the whole of the inner profile is affected by changes in 
suction velocity. This assumption will perhaps be quite acceptable when the inner 
layer is truly thin and when the suction varies slowly with position. In the present 
example, however, the application of a large amount of suction is probably respon- 
sible for the initial discrepancy, with the downstream improvement due to the 
outward diffusion effect. 


In spite of these criticisms it should be pointed out that the method predicts, 
with Head, that non-zero suction is required immediately, whereas the various one- 
parameter methods (e.g., Refs. 7, 8) would predict that zero suction is required 
initially. 


6. Applications to Some Flows with Zero Pressure Gradient 


In this section some examples are considered in which there is no pressure 
gradient, departures from the Blasius régime being determined entirely by the 
suction applied. In such circumstances, equation (23) for the skin friction may be 
simplified to give 


Suppose there is uniform suction downstream of a solid entry length /,; then 
equation (47) becomes 


T=1 when x<l.,, 


(1 32867 (T 1) T’ when x=> (48) 
a (49) 


then equation (48) becomes 


gv =0:32867 (7-1) (50) 


2 


This result indicates that the non-dimensional skin friction T is purely a function of 
t, =&'/7(1—£,/6)''*. If this were exactly true then it would be the same function of 
t, as in the case £,=0, studied by Iglisch°”. In fact it is not exactly true, but holds 
over a range of t, which increases with £,. 
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watson’) 
x PRESENT SOLUTION 


01 0-2 t 03 0-4 0:5 0-6 


Ficure 5. Skin friction for flat plate with uniform suction and solid entry length. 


Theoretically the range of validity of the present results is easily deduced from 
equation (27) and is 


Now Watson” has derived an analytic solution for this problem, which may be 
written in the form 


T = 1+ 2:21832, + 1:0215t,? — + (0:0446 — 0-3083£,-* *)t,*.. (53) 


This series supports the conclusion of the present method that, within the restricted 
range of values of f, given by condition (52), T is a function of t, only. Quantitatively 
it can easily be verified that the predictions of the present method agree well with 
the series (53) throughout the range over which each should be valid. As an 
example, the results for the case €,=1 are plotted in Fig. 5. We note that the 
agreement is within 6 per cent when t, < 0:4, and when f, =0°5 the agreement is to 
about 9 per cent. For values of f, greater than this value, Watson’s series does not 
appear to converge quickly enough to enable any firm conclusions to be reached. 


TABLE VI 
RESULTS FOR FLAT PLATE WITH UNIFORM SUCTION 
(SECTION 6) 


12 14 16 1:8 20 2°5 
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x PRESENT SOLUTION 


Ol o2 03 0-4 0-5 
Ficure 6. Skin friction for flat plate with uniform suction. 


Further, the results for the case €,=0 are shown in Table VI and compared 
with those obtained by Iglisch in Fig. 6. 
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A Theory for Flow with Axial Symmetry 
in Centrifugal Impellers’ 


R. J. TONKST, B.Sc., Ph.D. and A. G. SMITH, B.Sc., A.F.R.Ae.S. 


SuMMARY: A theory of incompressible flow with axial symmetry through the 
impeller of a centrifugal compressor is given more fully and more accurately 
that by previous writers. The theory is applied to compute the flow through 
two different impellers. For each design the experimental axial velocity distri- 
bution at the impeller eye is in good agreement with that computed numerically 
from the theory. This means that for the impeller assembly—not untypical 
of modern designs—considered here, theoretical calculation of intake axial 
velocity should provide a useful guide to designers. 


1. Introduction 


In certain respects one-dimensional theory does not describe the flow through 
a centrifugal impeller adequately. In particular, at the impeller eye, designers have 
been in the habit of modifying the blade angles predicted from one-dimensional flow 
theory, either to fit experimentally determined flow angles, or on a “cut-and-try” 
basis. Many writers have therefore felt it worthwhile to advance a little beyond 
one-dimensional flow theory. 


The actual flow in a centrifugal impeller is complex and, at present, intractable. 
It is usual therefore to make certain simplifications, namely, to assume the fluid to 
be inviscid and also to limit the flow to two dimensions. This two-dimensional 
limitation may be imposed in two ways : — 


(i) By assuming the flow to be completely constrained by the impeller vanes 
(i.e. axially symmetric flow) and then considering the flow referred to the 
meridional plane passing through the axis of the compressor (see Fig. 1(a)). 


(ii) By assuming the flow limited to two dimensions on surfaces of revolution 
generated by rotating the meridional streamlines of (i) about the axis of 
the compressor (the so-called “blade-to-blade” solution, see Figure 1(b)). 


*Most of this work was done in the course of Ph.D. research at City and Guilds College, 
London University. 


+Research Engineer, Austin Motor Company. 
Professor of Aircraft Propulsion, College of Aeronautics. 


Received April 1959. 
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Ficure 1. View of impeller channel showing : — 


(a) meridional plane, 
(b) surface of revolution, generated by rotating a meridional 
streamline about axis of impeller. 


These simplifications permit the description of the flow in terms of a stream 
function, the distribution of which can be determined with reasonable ease. From 
the distribution of the stream function, velocities may be deduced. 


A combination of (i) and (ii) constitutes the “quasi-three-dimensional” solution 
of the problem under consideration. 


Stanitz and Ellis“? have investigated theoretically the case of two-dimensional, 
non-viscous flow. The analysis was restricted to impellers in which the centre line 
of the passage between impeller blades generated a right circular cone when rotated 
about the axis of the compressor. The two-dimensional flow was then considered 
to lie on the surface of this cone. A further restriction confined the region investi- 
gated to that unaffected by the inlet configuration of the impeller and by the diffuser 
vanes. It is this latter restriction which is the main limitation on using the results 
for precise quantitative application. Later these two investigators extended their 
analysis to include channels whose centre lines generated arbitrary surfaces of 
revolution when rotated about the axis of the compressor. In this treatment a 
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Ficure 2. Cylindrical polar co-ordinates and velocity components relative to 
rotating impeller. 


knowledge of the streamline pattern in a meridional plane was assumed (i.e. an 
axial-symmetry solution had already been determined) and the mean streamlines 
were used to generate surfaces of revolution. Once again the effects of blade 
curvature at inlet were, of necessity, omitted. 


A rather interesting analysis developed by Stanitz and his colleagues“ 
considered axial-symmetry solutions for incompressible flow in an impeller fitted 
with curved inducer vanes. The analysis was limited to blade shapes which were 
composed wholly of radial fibres. With such blades no radial force is exerted on 
the fluid and the absence of this radial force simplified the analysis considerably. 


A comparison between a “quasi-three-dimensional” solution and a general 
three-dimensional solution has been made by the same author. He has concluded 
that the “quasi-three-dimensional” solution gives a sufficiently accurate description 
of the flow for engineering purposes. 


Marble and his co-workers“? have considered axially-symmetric flow through 
blades of arbitrary shape. Radial body forces arising from the blade geometry 
were taken into account. The results of his analysis are considered in Section 3. 


NOTATION 
r,9,z cylindrical polar co-ordinates 
R,9,Z_ dimensionless cylindrical polar co-ordinates (see Fig. 2) 
fi unit normal to vane surface 
absolute velocity 
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velocity relative to impeller 

angular velocity of impeller 

tangential velocity of impeller 

absolute vorticity 

vane (blade) force 

mean (dimensionless) axial velocity at intake 
fluid static pressure 
fluid mass density 
fluid total pressure ( 
p+ 4pV? 4pU" 


stream function, defined by equation (8) 


p+ 


dimensionless stream function 

stream function ratio vs 

angle defining the vane curvature (see Fig. 3) 
angle whose tangent is 7. / Me 

angle whose tangent is 7, / Me 

arbitrary points on vane surface 

relaxation net size 


residual errors associated with finite-difference equations 


components in the r-, 6-, z- directions, respectively 
impeller tip 

hub 

shroud 


grid points associated with relaxation net 


V.: VU, etc. 


2. Equation of Motion 


The equation of motion for incompressible flow, in vector notation, is 


February 1960 


- VH=v,w—v curl w+ F. 
ot op 


25 


U 
Ww 
F 
Vim 
H 
A 
P 
Suffixes 
|_| 


R. J. TONKS AND A. G. SMITH 


The equation of continuity is 


For steady inviscid flow with axial symmetry, the equation of motion reduces to 


In the present problem it is essential to work as many kinematic terms as possible 
into the equation of motion, since many of the boundary conditions are specifica- 
tions of fluid direction. For this reason the “vortex” form of the equation of motion 
(equation (3)) is the most useful, since it lends itself better to kinematic boundary 
conditions in the sense that v, w is essentially a velocity term and hence related 
directly to the boundary conditions. Further, (1/p)VH is often a simple quantity 
and indeed is zero in many problems. 


When rotary flows are considered, a more convenient, though less well-known, 
form is 


This has the virtue that in many rotary flows (1/p)V/ is zero and thus the 
equation only contains the terms F and Vw, both quite closely related to the 
kinematic boundary conditions. 


When 7» is zero upstream of an impeller, as it is for compressors without pre- 
whirl, V/ is also zero upstream. Since / is constant along a streamline in inviscid 
flow, if V/ is zero upstream it is zero everywhere, and equation (4) becomes 


Any method of throwing this equation into one in terms of the Stokes stream 
function must essentially incorporate relations between the tangential component of 
velocity and the axial and radial velocities, since the stream function itself is related 
to the axial and radial components only. These relations can be found by consider- 
ing that the relative velocity V is tangential to the surface of the vane; in fact 


V.n =0, 


whence (6) 

n,, Me and n, are geometrical properties of the impeller and may be obtained 
from the vane configuration (see Fig. 3). Use of equation (6) permits the left hand 
side of equation (5) to be thrown into “stream function” form. 
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It is as well to note at this juncture that 


F=n, 
so that F, = Fe = ; (7) 
n, ne n, 


1 ov 
V, and V,=- : 
ar r 0z (8) 
ov l 
ror r 
Now V.w.w=0. 
Therefore, from (5), F .w=0. 
Hence WM, + Welle + win, =0, 
n, n. 
Ne No 


Substitution from equations (6) and (8) into (9) yields 


~ 


Dy —20r™ ( ) ( Vi+r (V.)4 


Or? nme Cr Or\ ne none " Ne neecr 


n, 0 {n, n,n. Of n, 
+ = ( ) rv, = )+r (V,)+ r._(V.)+rV ) 
Ne cr\ Ne Me CZ \Ne Ne CZ Ne Ne CZ Ne CZ \M 


(10) 


In equation (10), which is the governing equation for the idealised flow, the 
first, third and fourth terms of the left hand side of the equation have been put in 
terms of the stream function, for reasons which will be apparent later. 


It is customary to specify blade shape in terms of certain angles rather than in 
terms of the components of the unit normal vector and it is necessary, therefore, to 
determine the relation between these components, which are needed for putting in 
equation (10), and the angles which will normally be known for a given design. 
These parameters are depicted in Fig. 3. The curvature of the guiding surface is, 
of course, purely arbitrary but, in the interests of clarity, has been depicted as being 
of constant chord, untwisted and extending to the impeller axis. This enables the 
angles A and » (essentially design parameters) to be clearly defined at any point P 
on the guiding surface. In fact such a surface is utilised in a common form of 
inducer vane, being very easily generated and consequently reducing the machining 
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MEAN SURFACE 


MEAN SURFACE 
OF BLADE OF BLADE 


my 
R 
/ AXIS 
/ 


ROTATION 


Ficure 3. Co-ordinate systems referred to mean surface of blade. 


difficulties in manufacture. For a twisted vane the angles A and v are not quite so 
readily accessible but, nevertheless, may be determined. It is easily shown that 


n, tana n, 
= and =ftanyv. . EE) 
Ne cosy No 


It will be seen that » defines the angle normally associated with radially- 


fibred vanes. 


3. Comparison with Existing Theories 


For impellers having radially-fibred vanes, n,=F,=0 and equation (10) 
reduces to 


(1 + tan? ») + 20r tan 
Substitution of tan tan «, 
Ne 


and use of the non-dimensional form, gives 


92%, (an? a 


+2Rtan<e=0, . (12) 
which is identical with the equation derived by Stanitz. 


Similarly, for axial turbo-machines—also for impellers fitted with inducer vanes 
sufficiently far upstream for them to be regarded as the initial row of rotating blades 
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At XX (upstream): — 
Distribution of is known. 
(ii) V, is constant over radial width. 


At YY (downstream) :— 
Distribution of is known. 
(ii) Vy, is constant over axial width. \ 


8 Cc 
INDUCER IMPELLER 


4 

foal 


In region A: 
ov 1 av 
OR? 02? R 
In region B: 
(1 + tan? ») + —— +tan v —_ (tan tan 0 
R Z eR 0102 V,,, 
In region C: 
l 
OR? eZ? R OR 


Ficure 4. Governing equations and regions of application 


in an axial compressor—the stream surfaces are practically circular cylinders and 
consequently the radial component of velocity may be ignored as a first approxi- 
mation in equation (10). Also, because of the small magnitude of the radial body 
force F,, it seems justifiable to approximate the local axial velocity V. (which is 
unknown) to the mean axial velocity V;,, (which will be a known design parameter). 


Thus equation (10) reduces to 


Oy day n, n, n, 0 ( n, ) 
=20r — + Vix 13 
which is identical with the equation derived by Marble’s method“. 


4. Simplifying Assumptions 


For the common type of impeller assembly used in the experimental verification 
of this theory, V, and v are both found to be small in the inducer vane region and 
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Ficure 5. Finite-difference approximations. 


hence it is reasonable to ignore terms in equation (10) which contain combinations 
of V, and tan v (or their derivatives) and also 0V,/0z, since these will be of third 
order of smallness. 


Thus equation (10) reduces to 


a7 


OY d*y 
(14) 


a 


(1 + tan? ») - —tan Ma 


r 


or, in non-dimensional form, 


(1 + tan’ + [ tan ap (lan w)—tany az (lan | aR az? —2R tan 
(15) 
Upstream and downstream of the impeller the tangential component of vorticity 
is zero and, therefore, 
1 ow 


It is convenient to introduce a stream-function ratio V=W/W, into these 
equations. The modified equations and their regions of application are shown in 
Fig. 4. 


5. Numerical Solution 


With the finite-difference approximations indicated in Fig. 5, equations (15) 
and (16) may be written, respectively, 


+ tan? IR, ~ 2 tan 37 (tan W, —2 (2+ tan?u,) V,+ 
h h 4R, tan po _ 
and ¥,+¥,+ (1- (1+ . (18) 
. 2R, 2 2R, 4 0 0 
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FiGure 6. Leading edge modification. 


It is a physical impossibility for a discontinuity to exist in the stream-function 
distribution at the leading edge of the inducer vane, although there is an abrupt 
change in the form of the governing equation at this point. Hence a method is 
required for blending the finite-difference equations (17) and (18) in the proximity 
of the leading edge, in order to ensure continuity of streamlines at this point. The 
method used was the derivation of a “weighed-mean” finite-difference equation 
applying at the leading edge. This is explained quite simply with reference to 
Fig. 6. 


Let F, and F, represent the residuals of equations (17) and (18) respectively. 


Let y be the significant fraction of F; applying at any given node; then (1 — y) 
is the significant fraction of F, applying at the same node. 


Hence the net residual at this node is 
F.= —y) 


The modified equation then becomes 


- h h 0 = 
[1 { tan’ 5 tan v, az (tan kz —2(2 + tan? 


[1+ 2R,*? { tan Hot tan vo 57 (tan Mo) 


(19) 
and this will apply on grid lines (0) and (1). 


Note: For the impeller design considered in this numerical solution. 
(R,? — Ry”)=0°102 
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For nodes lying along grid line (0), 


For nodes lying along grid line (1), y= ay 


Using equations (17), (18) and (19) a solution may readily be obtained, using 
relaxation methods, for the stream function distribution within an impeller, since 


the boundary conditions are known. 


6. Computation of Axial Velocity Distribution at Impeller Eye 


A numerical solution on the lines of Section 5 was obtained for a radial 
impeller with a tip radius of 4-00 in., rotating at 7,500 r.p.m. A mean axial velocity 
of 75 ft./sec. (nominal mass flow = 1,500 Ib./hr) was assumed upstream of the eye. 
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FIGURE 7. 


Inducer vane geometry. 
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~+--FLOW WITH ZERO TANGENTIAL VORTICITY 


—— FLOW WITH TYPE [ INDUCER VANES 


Radius R 


O95 OF O15 O2 OF O35 O4 OS OSS 
AXIAL DISTANCE Z 


Ficure 8. Streamline distribution in meridional plane—Type I inducer vanes. (Streamline 
designation indicates proportion of flow through impeller between streamline and hub.) 


Two types of inducer vane (see Fig. 7) were considered. Each was of constant 
chord and untwisted, the mean camber line being defined by 


J 


tan A= 


Both types were identical except that the leading edge of Type I was “cut back” 
in order to give approximate matching at the eye. It will be noted that such vanes 
are essentially non-radial and hence a radial body force is manifested. 


The fundamental flow between the hub and shroud of an impeller might be 
regarded as that with zero tangential vorticity everywhere. This is the flow that 
would be obtained if there were no blades on the impeller at all, and no pre-whirl 
in the air going through between hub and shroud. The distortion of the streamlines 
when an impeller and set of inducer vanes is present can then be judged relative to 
this fundamental case. A Zero-vorticity solution was therefore obtained in the 
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Ficure 9. Streamline distribution in meridional plane—Type II inducer vanes. (Streamline 
designation indicates proportion of flow through impeller between streamline and hub.) 


present case, since this provided a good basis for comparison with the non-radial 
inducer vanes. This solution is presented in Figs. 8 and 9, which show the stream- 
line distribution in the meridional plane. Similar solutions are presented for the 
Type I and Type II inducer vanes. 


Axial velocity distributions deduced from the stream-function distributions for 
Types I and II are presented in Figs. 10, 11, and 12. 


7. Experimental Measurement of Axial Velocity Distribution at 
Impeller Eye 


To confirm the theory developed, an experimental investigation of the air flow 
in the vicinity of the impeller eye was undertaken and compared with the theoretical 
conclusions. 


A 26-vaned radial impeller, fitted with 13 inducer vanes, was utilised and the 
two types of inducer vane discussed in Section 6 were considered. The compressor 
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was fitted with adjustable diffuser vanes which enabled the mass-flow rate, and 
hence the upstream axial velocity Vz,,, to be varied. 


The accurate measurement of the dynamic pressure so close to the rotating 
leading edges of the diffuser vanes presented certain difficulties which were overcome 
by the use of yaw-sensitive disc probes. The probes had previously been calibrated 
and were mounted on the intake shroud in such a way as to facilitate radial 
traversing and yaw determination. 


The rig layout is shown in Fig. 13 and the probe traversing mechanism can be 
seen in Fig. 14. 


Since the main preoccupation was the determination of the axial velocity 
distribution in the proximity of the eye due to the radial body forces imposed by the 
non-radial vanes, it was considered desirable to ensure that no velocity gradient 
existed by virtue of the upstream hub/shroud profile. A preliminary investigation 
was therefore undertaken far upstream from the eye; a radial traverse in the plane 
z= —0°875 in. showed the absence of any velocity gradient across the intake at this 
point (see Fig. 10). Similar radial traverses were made as close as possible to the 
leading edges of the inducer vanes (at z= —0-150 in.) for various values of upstream 
axial velocity. These demonstrated the presence of gradients for both types of 
inducer vane used (see Figs. 11 and 12). 


8. Analysis of Results 


From the results obtained in the investigation of flow conditions prevailing at 
the measuring station, situated 0-875 in. upstream of the eye, it was concluded that 
the assumptions of uniform axial velocity far upstream (see Fig. 10), also of zero pre- 
whirl, were both perfectly valid. 


The agreement between the mean axial velocity obtained from flow measure- 
ments at the Venturi meter and that obtained from measurements at the eye was 
quite reasonable (see Table I), there being less than 3 per cent difference in most 
cases. The very close agreement for upstream measurements (z= —0°875 in.) 
suggests that the difference is probably due to slight disturbances in the flow near 
the eye which were reflected in the probe readings. This is further verified by the 
fact that, for very low diffuser vane settings where there was a 6 per cent difference 
in the estimated values of mean velocity, wool-tuft observation indicated the 
presence of continuous blow-back at the eye tip. 


The absence of pre-rotation just upstream of the eye (due to drag from the 
rotating leading edges) was apparent from the negligible yaw angles obtained. 


The values of the velocity ratio included in Table II are also plotted on Figs. 11 
and 12. It can be seen that good agreement was obtained over three-quarters of 
the radial length for both types of inducer vane. Near the shroud, however, a 
disparity existed between the experimental and calculated values, the former being 
lower than anticipated. The experimental and calculated gradients are in good 
agreement, the plotted points lying above the calculated curves. This result was to 
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VARIABLE-DELIVERY COMPRESSOR 


Ficure 13. General layout of test rig. 


Ficure 14. Front view of compressor showing disc-probe transversing mechanism. 
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TABLE III 
DETAILS OF TESTS IN TABLES I AND II 


Test Air Flow (\b./hr.) Inducer Vane 


965 Type I 
1300 Type I 
1500 Type I 
1800 Type I 
1380 Type I 
1450 Type I 
1300 Type I 


be expected since the fall-off in velocity near the eye-tip results in a smaller mean 
axial velocity, and hence increased values of the velocity ratio Vz/Vz, resulting in 
the aforementioned displacement. 


Various values of mean axial velocity at intake were considered when using 
Type I inducer vanes. These ranged from 13 per cent below to 21 per cent above 
the value taken as a basis for the numerical solution. Such variations did not 
result in any substantial alterations in the experimentally determined gradients. This 
is useful, since it means that a numerical solution obtained for some basic through- 
put (provided that this is not too small) may be applied with a good degree of 
accuracy to other throughputs of about the same magnitude. 


The disparity between calculated and experimental velocities near the shroud 
was almost certainly due to disturbances in the flow originating within the impeller 
channels or even as far downstream as the diffuser vanes. Benser has confirmed 
that, at low capacity, flow reversal at the eye-tip originating inside the impeller 
channel may occur—a condition which was encountered by the authors when 
Operating the compressor at low diffuser vane settings and which represents the 
extreme case of flow disturbance inside the impeller influencing the upstream flow, 
particularly near the shroud. 


9. Conclusions 


The theory developed in Section 2 is considered verified in view of the results 
discussed. Although its development was restricted to idealised flow, it may be 
applied in the case of a real fluid flowing through an impeller in order, say, to predict 
the probable axial velocity distribution obtaining at the eye for a given set of design, 
and operating, parameters. In such a case, reasonable agreement may be expected 
over approximately three-quarters of the inducer vane height, but the velocity in the 
vicinity of the eye-tip is likely to be over-estimated and an inducer vane designed for 
the theoretical axial velocity will operate at a positive angle of attack. 


The numerical solution obtained demonstrates that distortion in the fluid stream 
will occur in the region of non-radial inducer vanes due to the imposition of radial 
body forces. Such distortion is unlikely to extend far beyond the inducer-vane 
region. Asa result of this distortion, a significant axial velocity gradient may occur 
at the eye, with consequent mis-matching of the fluid flow at the leading edges of the 
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inducer vanes, if these are designed on the assumption of zero axial velocity 
gradient. 
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Vibration of Rectangular Plates 
Point-Supported at the Corners 


HUGH L. COX and JACK BOXER 


(The Martin Company, Denver, Colorado) 


Summary : The fundamental frequencies of flexural vibration are determined for 
uniform isotropic rectangular plates that have free edges and pinpoint supports 
at the four corners. For the particular case of a square plate the lowest five 
frequencies and mode shapes have been determined. The second frequency 
is most unusual, since an infinite number of mode shapes exist with identical 
frequencies. Finite difference expressions, which simplify the treatment of the 
free boundaries for definite values of Poisson’s ratio, are used in conjunction 
with extrapolation procedures to obtain the approximate solutions. 


1. Introduction 


An exact solution of the differential equation of a vibrating uniform elastic 
plate that is point-supported at the four corners may be difficult to obtain; 
consequently an approximate solution has been determined through the use of finite 
difference expressions and an electronic digital computer. The finite difference 
approach was selected over other approximate methods such as Rayleigh-Ritz, 
collocation, influence coefficients, and so on, because the formulation time required 
for setting up the problem for machine solution is much shorter with the finite 
difference method. 


The usual assumptions regarding the bending of thin elastic plates are made. 
The effects of rotatory inertia and shear are assumed to be negligible. Poisson’s 
ratio » is taken as 0:3 for rectangular plates and 0 <= »< 0:3 for square plates. 
The plate under consideration is shown in Fig. 1. If the static lateral loads on the 
plate are replaced with equivalent inertia forces when the plate is vibrating freely, 
the governing differential equation for the plate is 


V2 V2w — w=0, y (1) 


where w is the plate deflection, D is the plate stiffness, p is the plate mass per unit 
volume, t is the piate thickness, and «, is a natural frequency in radians per second. 


Equation (1) may be expressed in finite difference form at a selected number of 
points. First-order operators, which include Poisson’s ratio, for points on and 
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FicureE 1. Types of plate considered. 


adjacent to the free boundaries have been tabulated by Jensen'’’ and Livesley and 
Birchall®. If equation (1) is written in difference form at each mesh point for a 
particular plate, the resulting simultaneous linear algebraic equations may be written 
matrically as 


where n is the plate length divided by the grid spacing (n=a/h), A is a matrix of 
coefficients of the V?V’w terms in equation (1), / is a unit matrix, W is a column 
matrix of deflections, and 


A= (3) 


Thus «,=n? J(A) =k, [(2,). 


Equation (2) is in standard eigenvalue form and the lowest roots A, may be deter- 
mined rapidly on an electronic digital computer for any representative matrix A of 
a particular plate grid. 


NOTATION 
Poisson’s ratio 
E Young’s modulus 
a,b plate length and plate width, respectively 
h_ grid spacing Ss Fig. 1) 
n=al/h 
w plate deflection 
D,p,t plate stiffness, density and thickness, respectively 


natural frequency (rad./sec.) 
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A, = ptw?a*/ D 
k, vibration coefficient, =n? /(A,) 
A,1I,W matrices defined just after equation (2) 


K elastic foundation modulus per length cubed. 


2. Fundamental Frequencies and Mode Shapes 


Plates with values of a/b of 1-0, 1-5, 2:0, and 3:0 have been considered. For 
each ratio of a/b the plate width has been taken as both four and six grid spacings, 
so that extrapolation could be used. Fig. 2 shows plan views of the plates that are 
six grid spacings wide. The fundamental mode shape will be symmetric about the 
two centre lines, as shown by the numbering system. Equation (2) has been solved 
on an electronic digital computer for each value of a/b, where b=4h and b=6h. 
Richardson’s h* extrapolation formulae“, 


36 (A,n*), — 16 (A,n"), 
36 — 16 = 
have been used to compute the vibration coefficients k, in equation (4). The 
calculations are tabulated in Table I. Values of k, for other values of nm lower than 
those given in Table I indicate that k, is approaching the true value monotonically 
in each case; consequently, extrapolation should give very accurate results. The 
values for k, given in Table I probably deviate from the “exact” answers by less 
than one per cent. 
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Ficure 2. Symmetric grid for fundamental modes. 
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TABLE I 
EXTRAPOLATION FOR FUNDAMENTAL FREQUENCIES 
a/b dint k, 
! 
4 0190281 48-7119 
6 0:038416 49-7868 50'6467 
6 0059568 77-1995 
9:69 92 
9 0011978 18:5868 190961 8-927 
8 0:020605 84-3988 
2 3 293 
12 0004123 85-4867 86°3570 
12 0:004225 87-6041 
q 
0000841 88-3230 88-8982 


The mode shapes for the grid system shown in Fig. 2 are given in Table II. 
Note that for the square plate the deflection at the mid-point of an edge is about 
two-thirds of the maximum centre deflection. For the other non-square plates it is 
interesting to note that the maximum deflection occurs at the mid-points of the long 
edges; a slightly arched condition occurs along the centre line parallel to the width. 


When a/b=00 the plate frequency reduces to the familiar simply-supported 
fundamental beam frequency for which k,=7*. Fig. 3 gives a plot of k, against a/b. 
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Ficure 3. Fundamental vibration coefficients for rectangular plates. 
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TABLE Il 
FUNDAMENTAL MODE SHAPES FOR PLATES SHOWN IN FIGURE 2 
(un =0°3) 


Grid a/b 
Location 1:0 1:5 2-0 30 


034407 0°35177 026141 0:17474 
2 0° 58825 065713 0°50327 0°34364 
3 067600 088141 0:70970 0°50168 
4 0°58026 1:0 086743 0°64423 
5 0°75981 0:11297 096631 0°76711 
082616 041240 1:0 0: 86668 
0°75981 067951 0:04188 094000 
089905 087885 028042 098489 
0:-95183 098513 0°50445 1:0 

10 0: 19186 069753 001135 
0:45973 0°84592 017900 
0:70293 093928 034210 
0°88653 0-97113 0°49540 
098502 0:07094 063408 
021992 0:29562 0:75387 
047731 0: 50849 0°85106 
071244 069315 0:92268 
089066 0°83569 096655 
098649 092561 098132 
095632 001919 

0:08124 0°18251 
0°30131 034197 

051040 04922 
069218 062855 
0-83274 0:°74642 
0:92148 084216 
0-95181 094275 
095600 
0:97057 
002196 
0: 18384 
0:34205 
049135 

0°62681 
074403 
083927 
090952 
095257 
0:96707 


3. The Square Plate 


A.similar approach has been used to find the first five frequencies and mode 
shapes of a square plate. Values of n=5 and n=6 have been taken so that extra- 
polation could be used. Fig. 4 shows the grid numbering system used for n=6. 
Table III gives the results of the solutions of equation (2) when »=0:3. Table IV 
gives the corresponding mode shapes for the grid shown in Fig. 4. The mode shapes 
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Ficure 4. Grid pattern for square plate for higher modes. 


have been plotted in Figs. 5 and 6 to illustrate the surfaces; dotted lines refer to 
node lines. 


It is most interesting and unusual that there are an infinite number of mode 
shapes with identical frequencies for the second frequency*. The lower part of 
Fig. 6 shows how two modes of the (a) type may be combined to obtain more basic 


TABLE Ill 
FREQUENCY EXTRAPOLATIONS POR A SQUARE PLATE 
(x1=0°3) 
Mode Extrapolation 
A k 
5 0362601 226625 
2 247-549 73 
6 0179798 233.019 
5 0453258 283.286 
325 
3 6 0:238087 308-562 19-13 
5 2°02722 1267-01 
6°40 4 
6 1-02700 133099 36°42 
5 239991 1499-95 
6 1-25090 1621°17 1996°68 od 


modes (b) and (c). The computer printed out the (a) type modes. As a check, 
case (b) was also obtained using finite difference expressions with a hinge (zero 
deflection) down the diagonal and case (c) was obtained with a hinge across 


*It could be possible that this phenomenon is a result of the numerical method of solution. 
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TABLE IV 
HIGHER MODE SHAPES FOR SQUARE PLATES 
(u=0°3) 
Grid Mode Number 
Location 2a 2b 3 4 5 

1 0°52207 042800 0: 50302 048159 0°82466 — 049310 
2 088521 0°79707 086736 85688 0-79791 -0°73475 

3 10 1:0 1:0 1:0 0 —0-78861 
4 084951 093765 0 86736 085688 —0:79791 — 073475 
5 0: 48397 0: 57804 050302 048159 — 082466 0 49310 
6 0:20274 — 42800 0:07502 —0-48159 082466 —0 49310 
7 047548 0 0°37920 0 10 — 035047 
8 066279 0 39652 060887 0°36362 0°73369 013858 
9 069381 069381 69381 0-50055 0 02813 
10 0°55495 0:82122 0 60887 036362 -0°73369 — 0: 13858 
il 0:28292 075841 0°37920 0 1-0 (35047 
12 — 005269 0:57804 007502 —0°48159 0 82466 0-49310 
13 029050 —(0°79707 0°07029 —0°85688 079791 73475 
14 0° 36694 0° 39652 021235 — (036342 0°73369 ~ 0 13858 
15 0 40496 0 032296 0 0 47978 0 44086 
16 0° 36438 036438 0: 36438 013585 0 069261 
17 0°24096 064592 032296 0 — (47978 0 44086 
18 0:05776 0:82122 0-21235 —0°36362 —0°73369 13858 
19 014994 093765 0.07029 —0 85688 79791 073475 
2 025390 10 0 10 0 --() 78861 
2 017616 069381 0 0: 50055 0 —0 02813 
22 009251 0: 36438 0 —0°13585 0 0 69261 
23 0 0 0 0 0 1:0 

24 009251 0°36438 0 — 013585 0 069261 
25 0-17616 069381 0 0° 50055 0 —0 02813 
26 025390 10 0 10 0 78861 
2 014994 —0: 93765 07029 085688 079791 —0°73475 
28 005776 082122 — 021235 0 36362 0:73369 — 013858 
29 024096 064592 — (32296 0 047978 044086 
30 0: 36438 0:36438 0: 36438 0-13585 0 069261 
31 040496 0 0°32296 0 047978 0 44086 
32 036694 039652 021235 —0°36362 0°73369 0 13858 
33 029050 0: 79707 0-07029 —0'85688 079791 
34 0:05269 0: 57804 0: 07502 48159 — 082466 49310 
35 — 028292 -0 75841 037920 0 —10 — 0-35047 
36 55495 —0°82122 — 0 60887 636362 ~-0°73369 0-13858 
37 — 069381 069381 —0-69381 0: 50055 0 002813 
38 —0°6627 039652 — 060887 036362 0°73369 0: 13858 
39 —0°47548 0 37920 0 10 —0°35047 
40) —0°20274 042800 —0-07502 — 048159 082466 -0 49310 
4l — 048397 —0°57804 — 0° 50302 048159 —0 82466 049310 
42 ~0°84951 093765 — 086736 0:85688 ~0°79791 —0°73475 
43 -1:0 1:0 0 ~0°78861 
44 —0 88521 -~ 079707 — 86736 085688 0:79791 —0°73475 


45 — 0° 52207 042800 —0 50302 0-48159 0°82466 — (49310 
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Ficure 5. Higher mode shapes. 


FIGURE 6. Some second mode shapes 
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FiGuRE 7. Square plate frequency variations with Poisson’s ratio. 


mid-span. These mode shapes agreed excellently (to six significant figures) with 
those obtained by the addition process of case (a) shown in the lower part of Fig. 6. 
Also, the frequencies agreed to eight significant figures. Of course, for the hinged 
cases (b) and (c) the first frequency corresponds to the second frequency for case (a), 
which is basically the problem under consideration. 


The first five natural frequencies for values of Poisson’s ratio of zero and 0-225 
also were obtained in order to construct Fig. 7. The fundamental frequency 
coefficient increases slightly with an increase in Poisson’s ratio; however, the other 
frequency coefficients decrease with an increase in Poisson’s ratio. 


Notice that the third mode shape of the plate under consideration is the same 
as the second mode shape of a free-free plate. Ritz“ obtained an upper bound 
k, =20-56 for this mode for a value of Poisson’s ratio of 0-225. The finite difference 
solution gave a value k, = 19-92 for this mode (see Fig. 7). There is approximately 
three per cent difference beween the solutions. This frequency decreases almost 
linearly with increasing values of Poisson’s ratio. This phenomenon has been 
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Ficure 8. Square plate frequency variations with Poisson’s ratio. 


observed in tests previously. The data shown in Fig. 7 is presented in a form 
compatible with most other published work on plates. It may be more reasonable, 
however, to factor out Poisson’s ratio from the plate stiffness parameter and to 
determine the frequency variation with Poisson’s ratio as (see Fig. 8) 


If the plate in Fig. 1 rests on a uniform elastic foundation, 


where K represents the elastic foundation modulus per length cubed. 
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The Pressure Distribution on a Wing with a 
Subsonic Leading Edge in the Presence of 
Anti-symmetric Body Shaping 


J. G. JONES 
(Royal Aircraft Establishment, Bedford) 


Summary: The pressure distribution on the wing of a centrally mounted wing- 
body combination due to anti-symmetrical body shaping (with respect to the 
wing plane) can be calculated in the case where the wing leading edge is super- 
sonic by using the fact that the upper and lower halves of the flow field are 
independent and that each half can be treated by quasi-cylinder theory. When 
the leading edge is subsonic the anti-symmetric body shaping produces flow 
round the leading edge. This paper is concerned with wing pressure distri- 
butions including the effect of this flow round the leading edge. 


1. Introduction 


The configuration considered is a lifting wing-body combination comprising a 
wing, with a leading edge subsonic over at least part of the span, centrally mounted 
on a quasi-cylindrical body. The body shape is obtained by adding small distor- 
tions to a basic circular cylinder; the distortions may be considered as the sum of 
two parts, one symmetrical and the other anti-symmetrical with respect to the 
wing plane. 


This paper is concerned with the pressure distribution on the wing due to the 
anti-symmetric part of the body distortion, including the interaction with the 
subsonic leading edge. 


Within the limitations of linearised theory the effects of anti-symmetric body 
distortions are independent of the effects due to 


(a) wing thickness, camber, and incidence, 
(b) body incidence, 
(c) symmetric body distortions. 


The pressure-fields due to all these effects are simply additive. Consequently 
it is sufficient to consider a quasi-cylindrical body with anti-symmetric body 
distortions in combination with a flat plate wing at zero incidence. The results are 
then applicable to more general wing-body combinations. 
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Figure 1. General configuration. 


The configuration is illustrated in Fig. 1. In the first place it is assumed that 
the upstream limit AB of the pressure field due to body shaping meets the leading 


edge at a point B away from the body side. This implies that either the leading 
edge is supersonic at the root or that the body shaping starts downstream of the 
root leading edge. Consideration is later given to an approximate solution when 
this condition is violated. 


Upstream of BC the pressure on the wing is that due to the direct effect of body 
shaping on the same side only. This pressure field is given by quasi-cylinder 
theory”. 


Downstream of BC the pressure on the wing is influenced by the body shaping 
on both sides of the combination. 


Downstream of BC, but upstream of CD, the pressure on the wing is that due 
to the direct effect of body shaping on the same side, together with the effect of the 
upwash field between the leading edge and the upstream limit of the effect of body 
shaping. The theory of this paper applies exactly to this region and involves the 
application of Evvard’s method? to the upwash field. (The treatment of Evvard’s 
method is derived from that in Ref. 3.) 


The solution of problems of this type has been generally discussed in Ref. 4. 
This paper is concerned with two particular cases in which explicit expressions 
involving tabulated functions for the pressure distribution due to the upwash field 
can be obtained. The two cases are 


(i) body distortions in which the cross sections are circular, 
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FiGure 2. Illustration of theory. 


(ii) body distortions of arbitrary cross section, but with the restriction that the 
leading edge is “nearly-sonic” (this expression is defined precisely later). 


Downstream of CD the reflection by the body of the pressure field due to the 
upwash must be taken into account. The problem of finding the reflected pressure 
field is not treated here. However, an estimate of the order of magnitude of this 
reflected pressure field can be obtained by finding the flow due to the upwash across 
the body side downstream of C in the wing plane. 


The effect of an upwash field on the opposite wing involves diffraction of its 
pressure field around the body. This problem is not treated here. In many 
practical problems the effect is negligible. 


A particular example of the case in which the subsonic leading edge extends 
right into the body side and anti-symmetric shaping starts at the root leading edge is 
discussed in Section 3. The leading edge is straight and swept behind the Mach 
line by a small angle, and the pressure coefficient on the wing is found to the first 
crder in this small angle. Conditions at the wing root leading edge are discussed 
in detail. 


NOTATION 


}rectangula co-ordinates 


6 meridian angle 
£,» characteristic co-ordinates 
¢=£/n, conical variable 
@ disturbance velocity potential 


?, (x, y)=¥, (€,), velocity potential, deflected wall 
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j 
! 
as ) 
fet 
2€.7,) 
«3 
(6,7), 
4 
(x,7,:; 


¢, (x, r)=¥, velocity potential, quasi-cylinder theory 
9. (x, r)=¥, velocity potential, flow round leading edge 
t(x,r)=2(€, »), upwash distribution 


n=F (£), equation of leading edge 
A.,(x) multipole strength 
H(y) function illustrated in Fig. 2 
M.,,(x) influence function 
R_ body radius 
7 parameter in definition of body shaping 
2.n(x) function in definition of body shaping 
free stream velocity 
W..(x,r) influence function 
P= —(2/V) 0?/0x, linearised pressure coefficient 
P, pressure coefficient, quasi-cylinder theory 
P, pressure coefficient, flow round leading edge 
a;,b; functions in definition of pressure coefficient 
a angle defining leading edge sweep 
S,5,,5, regions of upwash field 


=,=,,2, regions of upwash field in characteristic co-ordinates 


2. General Theory 
2.1. THE APPLICATION OF EVVARD’S METHOD 


Figure 1 illustrates the configuration considered and the cylindrical co-ordinate 
system. Anti-symmetric shaping starts at x=0. For simplicity, restriction is made 
to the case 8=(M?-—1)'/?=1 throughout. Other Mach numbers can be treated by 
the usual transformation known as Gothert’s rule. 

At present attention is confined to the region BCD. Then the disturbance 
velocity potential <P (with u=V®) may be expressed as the sum of the potential ¢, 
due to the shaping on the same side of the wing, given by quasi-cylinder theory, and 
the potential ¢, due to a source distribution over the upwash field, given by 
Evvard’s method. 


Then ©® is anti-symmetric and continuous upstream of the leading edge. 
So, over the upwash field, ‘ (1) 


Characteristic co-ordinates are now defined by 


(2) 
y=xXt+re+1 


The body radius is taken as the unit of length. 
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Consider a point P (x,, r,; €,, »,) on the wing and suppose 


(x,, r,)= v m) 
(,, 


The upwash distribution is denoted by 


t (Xo, =8 (Eos No)» 
Then ®(x,,7,)=9%, 1.7), ©) 
(S is the part of the upwash field within the upstream Mach cone from P.) 
In terms of characteristic co-ordinates (6) becomes 
1 dé, { (Eon) dtm} 


= Fé.) 


where »=F(€) is the equation of the leading edge and the other quantities are 
defined in Fig. 2. 


We now consider a point Q in the upwash field on the Mach line y=», through 
P. Qhas characteristic co-ordinates »,) where O< & < €”. 


From equation (1), we have 


Po 7,)=0. 


Then relation (5) gives 


4 


(n, 
) 


for allO< < &”. 


Equation (8) may be regarded as an Abel integral equation of which the 
solution is 


(0, No) dno d ¥, (t, dt 9 
ag, | (§ ” 
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Substitution of equation (9) in equation (7) gives ¢, (x,,r,)(=v.(£,.,)) in 
terms of integrals over the upwash field. 


The rest of this paper is concerned with cases in which the expression for 
$, (x,,7,) takes a particularly simple form. 


2.2. REPRESENTATION OF BODY SHAPING BY SOURCES AND MULTIPOLES ALONG 
THE Bopy 


Quasi-cylindrical body shaping can be represented by sources and multipoles 
distributed along the body axis“. We consider the representation of the shaping 
on the upper half of the configuration. Body shaping, and corresponding pressure 
distributions, on the lower half are equal in magnitude and opposite in sign. 


Then it is sufficient® to consider only multipoles of which the potential is 
given by 
z,-f, 


-1 
where @ is the meridian angle and equals zero in the plane of the wing, and 
Aon ( = 1)=0. 
2.3. BODY SHAPING WITH SEMI-CIRCULAR CROSS SECTIONS 


In this case the body shaping can be represented by a line of simple sources. 
The potential is then given by 


1 [  A,(x,) dx, 
In terms of characteristic co-ordinates 
where €,=x,+1. 
Then equation (8) becomes 
" 


0 
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Equation (13) may be regarded as an Abel integral equation of which the 
solution is 


(€,) g No) dn, 


=0 
(n, £,)*/? (n, 
F¢é,) 


for O< 
Substituting from equation (14) into equation (7) 


$2 
H(y,)-1 
l j A, (x,) dx, 


-1 


where the function H (n) depends on the shape of the leading edge and is illustrated 
in Fig. 2. 


Substitution from equations (11) and (15) into relation (5) gives 


__A, (%,) dx, 


® (x,, | {(x, —r,?}} 2° . (16) 
H(y,)-1 


Equation (16) expresses a result exactly analogous to Evvard’s result for a 
planar distribution of sources. In words: “the potential at P due to the upwash 
field is cancelled by the potential of the linear distribution of sources that produces 
the upwash field.” The potential at P is just that due to the remaining source 
distribution. 


This result does not hold for a distribution of higher multipoles. However, 
for a combination with a “nearly-sonic” leading edge, simple analytical results are 
still obtained (Section 2.6). 


2.4. DISTURBANCE VELOCITIES 


In this section are presented the values of 09,/0x, (from which the linearised 
pressure coefficient due to the upwash field is deduced) and [¢¢,/¢r,] Ni , (which 


gives the flow across the body side due to the upwash field). 
Equation (15) gives, after simplification, 


H(n,)-—1 
09, SLA, = A, (H (n,)- 1\(1 — H’ (,)) AS dx, 
- 
{(x, — x.) } 
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or, = 


{(x, — (n,)+ 1 = 
H(y,)-1 
\ Ag’ (x(x, dx. 


where, in equation (18), 7,=x,+2. (It is convenient to integrate equation (15) by 
parts before differentiating.) 


Equations (17) and (18) hold for H(n,)>0. For H(n,)<0, both 09,/0x, 
and vanish. 
Note that, for points (x,, r,) on the leading edge, £,=x,-—r,+1=A (n,). 


So 09,/0x,, and hence the linearised pressure coefficient, has a square root 
singularity on the leading edge for H (n,) > 0. 


2.5. RELATION BETWEEN SOURCE STRENGTH AND BopDy SHAPE 


The relation between source strength and quasi-cylindrical body shape with 
semi-circular sections is given in Ref. 1 and is 


z, 


Ay | 9, (0M, dt, 


0 
where the body shape is given by 


OR 


and M, (t) is an influence function, with a square root singularity at t=0, which is 
defined and tabulated in Ref. 1. 


The most general form of g,(x) to be considered here is a continuous and 
differentiable function, except possibly for a finite number of discontinuities. 


Then, from equation (19), 


z, 


A, (x, —1) dg, (1), ‘ (21) 


- 
where the integral is taken in the Stieltjes sense. 


If g, (x) has its only discontinuity at x=0, equation (21) can be written 


Ac! (x) () My at. (22) 
0 
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Equations (19) and (21) can now be used with equations (17) and (18) to give 
the disturbance velocities due to the upwash field in terms of the body shape. 


In the special case where 


x)=0, x<0 
(23) 


=g, (constant), x>0 


(i.e. g, (x) is a multiple of the Heaviside unit function, and the body shaping consists 
of a straight taper starting at x =0), we have, from (22), 


A, (x, —1)=22V7M, (x2) 
and equations (17) and (18) become 


H(»,) 


Ox, i(x, —H (y,)+ 1)?—r,?7}? 2 M, dx, 


H(»,) 

_M,(x,) dx, ] 
(x, —x, + 
0 


M, (x.) dx, + 


[ ar, por {(x, —H (y,)+ 1)? 


H(, 


| (x, —x,+1) M, (x,) 


where, in equation (25), »,=x,+2. Relation (5) gives the full expression for the 
disturbance velocities, 


OP (x,,7,) 06,(x,.7,) 00. (%,.97,) 
Ox, Ox, Ox, 


with a similar expression for [0 (x,,r,)/0r,],,<1. 


This section is concluded by quoting the results“ of quasi-cylinder theory for 


substitution in relation (26): — 


~ 


Ox, 
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Ficure 3. First illustrative example. 


Here £,=x,-—r,+1 and W, (t,r) is an influence function defined and tabulated in 
Ref. 5. 


The results of the application of the foregoing theory to the example illustrated 
in Fig. 3 are shown in Fig. 4. For this example H (y)=(y-—4)/3. The anti- 
symmetrical shaping is taken to be defined by conditions (23) on the upper surface. 
The results are expressed in terms of linearised pressure coefficient as follows : — 


2 OP (x7) _ 


P= - 
V Ox, : (29) 
P,=—- 2 09, r,) P,=- 093 


2.6. “‘NEARLY-SONIC” LEADING EDGE 


The second case in which relatively simple analytical expressions for the 
disturbance velocities due to the upwash field can be obtained is for combinations 
with “nearly-sonic” leading edges. This concept is illustrated in Fig. 5. For all 
points R on the leading edge downstream of the line AB, 


tany <1. ‘ . (30) 
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It then follows that for any point (x,,r,) in the upwash “eld 


x,+1-r, 


We now consider anti-symmetric body shaping which is given by 


OR 
- = < 6 
(x) cos 2n 


on the side on which the velocities are to be found. 


2dr/ox 3 DR/d 2 
0-8} Qer 
0-6} 
0-4) 
O2}+ 0-2} 
f 
LEADING MACH MACH MACH CH magn 
EDGE LINEAB LINE BC LINE CD LINE AB LINE CO 
AND LEADING EDGE 
2oR/ox ul 
(a) r<3 
re 
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> 
/ 
oO 
/ 
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FiGURE 4. Wing pressure distributions. 
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Ficure 5. ‘“Nearly-sonic” leading edge. 


At a general point (x,,7,) of the upwash field the potential due to body 
shaping is 


_ cosh [2n cosh-* {(x, x2)/ro}] dx, 


where Aan (X-1)=28Vr San (t) Man —1) dt 


and M.,, (t) is an influence function defined and tabulated up to n=3 in Ref. 1. 


Under condition (30), using condition (31), 


x,t+1l—r 


0 


and so, to a first approximation, 


cosh 2n cosh~? (* = cosh [ 2n cosh-! (1 ~ 1. 


A first approximation to equation (32) can therefore be taken as 


N 
1 Aa ax. A: (x,) dx, 
9, 
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Equation (34) shows that over the upwash field the potential due to the multi- 


pole distribution can be approximated by that due to a simple source distribution of 
strength = A,, (x). 
n=0 


In this case equation (8) becomes (cf. equation (13)) 


0 F ¢,) 


for allO=<=& =<”. 


Equation (35) may be regarded as an Abel integral equation of which the 
solution is 


N 

= é " 
+ 8 (£0.70) _ (36) 
(n, —€,)'? (n, —1,)'!? 


for all O< &”. 


Substitution from equation (36) into equation (7) gives 


From equation (37), equations can be obtained for the disturbance velocities in 
exactly the same way that the equations (17) and (18) were obtained from 
equation (15). Equation (33) can then be used to express the results in terms of 
body shape, exactly as relation (19) was used in Section 2.5. 


3. Approximate Solution when the Upwash Field Extends to the 
Body Side 


3.1. GENERAL REMARKS 


Referring to the general arrangement (Fig. 1), it was pointed out in the 
Introduction that the method given in this paper only applies exactly to points on 
the wing upstream of CD. When the upwash field extends right into the body side 
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Ficure 6. Second illustrative example. 


the point C coincides with A at the wing root leading edge and the region of exact 
application of the theory disappears. 


This case will be discussed by means of a particular example, illustrated in 


Fig. 6. The wing leading edge is straight and swept behind the Mach line by a 
small angle a (a particular case of a “nearly-sonic” edge). For this example 
H (»)=(n—2)tan 2. The anti-symmetric shaping is taken to be defined by (23) on 
the upper surface. 


The linearised pressure coefficient at a point on the wing of this configuration 

is now expressed in the form of an infinite series 
P=P,(l+a,+a,+.. .)+P,(1+6,+6,+...), (38) 

where the a,, b;, as well as P,, P,, are functions of position and where :— 

P, is given by quasi-cylinder theory, neglecting the flow round the leading edge, 
P, is given by the method of this paper and is due to the flow round the leading 
edge associated with P, but neglecting reflection by the body, a,P, is given by quasi- 
cylinder theory and is obtained by cancelling the flow across the body cylinder due 
to the flow field associated with P,, b,P, is obtained from a,P, in the same way that 
P, is obtained from P,, and so on. P, has a square root singularity along the wing 
leading edge. 


It is shown in Section 3.3 that relation (38) can be written 
P=(P,+P,X1 +O (2*"*)). (39) 
Neglecting terms of O (2"'*), we obtain the approximate solution 


P= P,+P,. . ‘ ‘ (40) 
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THE 2-AX'S NORMAL 
TO THE PLANE OF THE 
PAPER 


Ficure &. Illustration of solution at root leading edge. 


This approximate equation holds all over the wing. Near the wing leading 
edge P, is the dominant term. However, in the wing root P, is small, in fact of the 


same order as P,a,=O(2°'*) (since P,a, is obtained by the reflection of P, by the 
body), and so in the wing root equation (40) is replaced by 


The approximate solutions (40) and (41) are illustrated for several stations on 
the wing in Fig. 7, for the case tan 2=1/6. 


At the wing root leading edge conditions become conical and it becomes 
possible to find higher order terms in the series (38). This is done in Section 3.2. 


The results of Section 3.2 are used in Section 3.3 to estimate the errors in the 
approximate equations (40) and (41). 


3.2. THE SOLUTION AT THE WING Root LEADING EDGE 


This section is concerned with the conical configuration illustrated in Fig. 8. 
The origin O is chosen at the wing apex. For y > 0 the configuration corresponds 
to conditions very near to the apex A in Fig. 6. The curved basic body side of 
Fig. 6 is replaced by a plane wall, y=0. The body shaping, 0R/@x=7g,, for x > 0, 
is replaced by a constant deflection of the wall, éy/éx=rg,, for x >0,z>0. For 
z <0 the wall deflection is equal in magnitude and opposite in sign. The linearised 
boundary condition is applied in the plane y=0. 
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For z>0 the boundary conditions are satisfied by adding to the (two- 
dimensional) potential due to the deflected wall (which would be the complete 
disturbance potential if the wing leading edge were supersonic) a distribution of 
sources over the region S,, in the plane z=0, between the leading edge and the 
leading Mach line, and an (equal) image distribution for y < 0 over the region S, to 
satisfy the boundary condition on y=0. 


The boundary condition to be satisfied on the upwash field is 


where ® is the disturbance velocity potential. 
This leads to a problem of a type discussed by Evvard™, involving interactirz 
upwash fields.* 


In terms of the Cartesian co-ordinate system illustrated in Fig. 8, the following 
characteristic co-ordinates are chosen :— 


| 
t 


The (two-dimensional) potential due to the deflected wall is denoted by ®,. 


Suppose the upwash distributions over S,, S, (~,, =, in characteristic 
co-ordinates) are 


t, (x, 9) | 


t (x, (En) J 


2, (é,)=2, ©). . (44) 
Then, at P(x,, 
2m {(€, — — {(€, — — ) 


2 


*The problem, being conical, could no doubt be solved using complex variable techniques. 
The solution given here is more suitable for the present purpose of discussing the terms in 
equation (38). 
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In particular, 


7) z 


1 
dé, | g dn, 
) — 
 (x,, 0) (x,, 0) + (x, . (46) 


te 


At a point (£,, ,) in the upwash field S,, we have, combining equations (42) 


and (45), 
‘a 
¢ 82 No) "o (Cas? 


for all O< €, Sy, tana. 


Equation (47) may be solved as an Abel integral equation, using relations (43) 
and (44), giving 


fy tana 
81 (nos &) dn, 2, No) dn, = 


cot a 


4Vg,7é,''?. . (48) 


Equation (48) may itself be regarded as an Abel integral equation, thus 
reducing to 


1 0 | dn, { 1/3 (7, 
cota 0 


0 


fy tana 
cota2)'/? «=(n,—&, cot 2)'/? 


Introducing conical variables 
and writing g, (é,.,)=2, we have, for Stan 2, 
tare 


4V 2,75)? { (p) (cot — p)' (49) 


8: = =f, cot AV gor 


Equation (49) is solved by successive substitutions, writing 


4V 2,74," 
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+ 


x J(1—pcot a)'/? (1 — pt.) 


0 


where ((,)= — 


(51) 


taba 


It can be shown by induction that 


an2 n 
lem < tan?’a 
and so, for z small enough, the solution (50) converges faster than some 
geometric series. 


Substituting from equation (48) into equation (46), and expressing the result 
in conical variables, we obtain 


tang lang 


x dp . g, (s) ds 
(x,,0) 1 (%,,0)+ = J(1—p)'!? P (1 — ps)'!? (92) 
0 0 
Substitution of the solution for g,(¢,) given by equations (50) and (51) 
gives finally 


OP (x,,0) (x,, 0) 
CX, Ox, 


(tan + (tan — (tan 2)'/? +O (tan ayia} 
(53) 


and hence 


2 oP (x,, 0) 22,7 { \ 
V ae. 28 3% O (2°!) (54) 


P= 


This, then, is the pressure coefficient at the wing root leading edge of the 
configuration in Fig. 6. 


The successive substitutions correspond physically to successive reflections in 
the plane wall, and so the result, equation (54), can be interpreted in terms of 
equation (38). At the wing root leading edge, where the body reflects as if a plane 
wall, P.=P,a,, P.b,=P,a,, . . . and soon. Then the first term in (54) is P,, the 
term in 2°’? is the first term in the expansion of P,+P,a,, and higher order terms in 
(54) successively get contributions from P,b,+P,a,, . . . and so on. 


3.3. ERRORS DUE TO APPROXIMATION 


This section is concerned with the errors in the approximations (40) and (41) 
to equation (38). 
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To find a,P, (equation (38)) exactly it would be necessary to find the flow across 

the body cylinder due to the flow field associated with P,, which is represented by a 

source distribution over the upwash field. An estimate of the magnitude of this 

cross flow can be obtained by finding the flow across the body cylinder in the wing 

plane only. This is given vy equation (25). Using the series expansion of M, (t) 
given in Ref. 1 it can be shown that, for small 2, equation (25) is equivalent to 
2/2 


To the first order in 2 this function is greatest at x,=0 and decreases mono- 
tonically as x, increases. At x,=10 it has fallen to 0:32 of its value at x, =0. 


Taking this as an indication of the general behaviour of the flow across the body 
cylinder it follows that a, is O (2°/*) and decreases (for a fixed value of r,) when 
x, becomes large. 


b,P, is obtained from a,P, in the same way that P, is obtained from P,. So b, 
is also O (2*/*) and decreases for large x,. 


In the wing root equation (40) is replaced by equation (41). The neglected 
term P, decreases with increasing x, (Fig. 7). So the error in approximation (41) 
decreases for large x,. Use of equations (53) and (54) gives the error in the 
approximation (41) at the wing root leading edge. For example, when tan z=1/6, 
the error is approximately 6 per cent. This will be typical of the error in the wing 
root for small x,. 
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Struts with Linearly Varying Axial Loading 


J. S. PRZEMIENIECKI, Ph.D., A.F.R.Ae.S. 
(Bristol Aircraft Limited) 


SUMMARY: Analytical solutions for the deflections and bending moments are 
derived for elastic struts with linearly varying axial loading applied either 
(a) along the undeformed straight strut axis or (4) at the strut axis in the 
deflected position. The instability of struts with various end conditions is 
discussed and the stability criteria are given as a series of curves relating the 
maximum compressive axial load with the maximum distributed loading 
Furthermore, explicit formulae for the deflections and bending moment distri- 
butions have been compiled for several typical cases of lateral loading 
combined with the axial loading of the type (a), which is of some practical 
importance in aircraft structures. A numerical example is included to show 
the practical application of the general analysis to a typical strut problem. 


1. Introduction 


Struts with linearly varying axial load can be divided, in general, into two basic 
types, depending on the method of application of the distributed axial loading, 
which may be applied either (a) along the undeformed straight strut axis or (5) at the 
strut axis in the deflected position (see Fig. 1). 


The strut problems in aircraft structures can often be approximated to the 
type (a); the distributed axial loading can be taken as a constant shear flow q acting 
in the plane containing the undeformed strut axis, but once any flexural deforma- 
tions have taken place this loading becomes equivalent to a distributed loading q 
applied at the strut axis in the deflected position and a distributed couple n=qy, 
where y is the flexural deflection of the strut. The type (b) occurs in struts loaded 
by their own weight, where the distributed loading q is applied at the strut axis in 
the deflected position. 


The analysis of struts with uniformly distributed loading applied along the 
undeformed strut axis does not appear to have been published previously. On the 
other hand, struts with distributed loading applied at the strut axis in the deflected 
position have been analysed quite extensively in technical literature; the earliest 
paper on this subject is due to Greenhill, while a more general theory has been 
developed recently by Biiltmann®. 


The main object of this paper is to present a general analysis of uniform elastic 
struts subjected to linearly varying axial loading, with particular reference to the 
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TYPE (a) TyPE (b) 


FiGureE 1. Basic types of axial loading. 


two types of loading discussed earlier. The mathematical formulation of the 
problem and its solution present no difficulty. The derivation of the governing 
differential equations for the flexural deflections of struts under a general loading 
has already been given by Duncan™. The differential equations for the deflections 
can be solved with the aid of Bessel functions of order 1/3 but the practical 
application of these solutions is somewhat unwieldy. Another alternative is to use 
solutions in series form which can then be used directly for the calculation of 
deflections, slopes and bending moments. Such solutions would also allow a direct 
comparison with the special case of constant axial load which otherwise would not 
be possible with the solutions in terms of Bessel functions. In the present paper the 
series solution has been used as a basic solution, but the actual computation of the 
deflection functions and their first derivatives has been carried out via the Bessel 
function solution, because of slow convergence of the series solution; thus the 
necessity of computing a large number of terms in the series solution has been 
obviated. Only for the special case of a strut with zero axial load at one end, has 
the Bessel function solution been used directly. 


The stability criteria for struts with various end conditions are given as a series 
of curves relating the maximum compressive axial load with the maximum 
distributed loading. This shows that struts with distributed loading applied along 
the undeformed strut axis are more stable than struts with loading applied at the 
strut axis in the deflected position. 


Furthermore, explicit formulae for deflections and bending moment distribu- 
tions have been compiled for several typical cases of lateral loading combined with 
axial loading of the type (a), which is of some practical importance in aircraft 
structures. 
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A numerical example is also included to show the practical application of the 
present analysis to a typical strut problem. 


NOTATION 


y, flexural deflections 


x distance measured along the loading axis 


1 strut length 
e eccentricity of loading 
I moment of inertia of the strut cross section 
E Young’s modulus 
M_ bending moment 
.M, bending moments at x=0 and x=/ respectively 
S shear force 
Q lateral concentrated load 
P=P.,+ qx, axial loading 
P, axial load at x=0 
q_ distributed compressive loading 
w distributed lateral loading 
n distributed bending couple 
n. constant distributed couple 
=P. / El 
z= {2/(32°)} {1+ 
Ze = 2/ (30°) 
6=(2/3Ax)? 
Yer complementary functions 
..¥, complementary functions of equation (8a) 
"1, %, Complementary functions of equation (17) 
y, particular integral of equation (8a) 
», particular integral of equation (17) 
yp, integral defined by equation (57) 


®, (ux), P, (ux) Bessel function solutions of equation (44), defined by equations 
(46) and (47) respectively 
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C,,...C, constants of integration 
A,...,D constants of integration 
a,b,c,d constant coefficients 
k,,k, constants 
W =P /P, Wronskian of equation (44) 
J.( ) Bessel function of the first kind of order » 
Y.( ) Bessel function of the second kind of order »v 
) gamma function 
F (ux) function of ux 
o=dy/d (ux) 
f(x), ¢(x) functions of x. 


2. Governing Differential Equations for Defiections® 


Consider the equilibrium of the element AB of length ds in a uniform strut 
subjected to distributed loading and let Ox be the axis of the strut when it is 
unloaded (see Fig. 2). The distributed loading consists of lateral loading w, 
compressive loading g and bending couple n; both w and q are assumed here to be 
constant. In the notation used in Fig. 2, P denotes the total load in the direction of 
xO applied to the right of the point A, while S is the total load in the direction of 
Oy applied to the right of A. The bending moment M at any cross section of the 
strut is the clockwise moment due to all external loads between x and x=/, where / 
is the strut length. 


Taking moments about the point B and neglecting second-order terms, 


dy 


Figure 2. Forces and moments on element AB; for loading type (a), n=n,+qy, while, for 
type (6), n=n,, where n, is a constant distributed couple. 
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Resolving forces along the y and x axes respectively, 


dS 
dx 


dP 


and 


Equation (3) becomes, on integration, 
P=P,+ 


where P, is the axial load at x=0. 


From the linearised simple bending theory 


M=—E1 22. 
ax? 


(5) 
where E/ is the flexural rigidity of the strut. Differentiating equation (1) with 
respect to x, and then substituting equations (2), (4) and (5), the governing 
differential equation for the deflections is obtained: 


de { EI de (P,+qx)y ft +4 dx 


For q=0 and n=0, equation (6) reduces to the well-known strut equation 


dty 


El w. ‘ (7) 


Solutions of equation (6) are known for n=0, i.e. for struts loaded by distributed 
compressive axial loading applied at the strut axis in the deflected position. An 
important case, however, arises in the analysis of compression members attached to 
shear panels. This may be represented approximately by a strut with linearly 
varying axial loading equilibrated by a distributed shear (shear flow) acting along 
the undeformed straight strut axis. This distributed shear, when transferred to the 
neutral axis of the deflected strut, becomes equivalent to a uniform distributed 
loading and a distributed couple n= qy. 


A constant distributed couple m. can also be conveniently included in the 
general analysis, so that for the two basic types of loading the distributed couples 
can be iaken as 


(a) n=n.+qy and (b) n=n.. 


In Sections 3, 4 and 5, a detailed discussion of the general solution of 
equation (6) is given, with particular reference to these two types of loading. 
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Substituting n=n,+qy in equation (6), the following basic equation is obtained 


d? 


dx? 


{EI dx (P,+ qx) y we. (8) 


Equation (8) does not involve ,., but when this equation is solved n, will appear 


in the constants of integration calculated by satisfying the end conditions at x=0 
and 


Introducing now two non-dimensional parameters »?=P,/EI and z=A/n, 
where A° = q/ EI, equation (8) can be re-written as 


d (ux)? +(1+ a°ux) y 


Integrating this equation twice, it follows that 


w 


y" + (1+ px) y=Cy+Cyix+ 5 "op 


er. . (9) 


where C, and C, are arbitrary constants of integration and primes (’) denote 
derivatives with respect to px. 


The complementary functions of equation (9) can be obtained in a series form, 
using standard methods, so that 


where C, and C, are arbitrary constants, while y, and y, are the two complementary 
functions. The detailed analysis used to obtain the solutions y, and y, is given 
separately in Appendix I, where the series solutions are derived in terms of Bessel 
functions. 


The particular integrals of equation (9), together with the solutions y, and y,, 
give the four complementary functions and one particular integral of equation (8a). 
Thus the complete solution of equation (9) or (8a) can be written as 


where y,, y, and y, are the particular integrals of equation (9) with the right hand 
side equal to 1, (ux) and (x) respectively. 

The solutions for the particular integrals of equation (9) are derived in 
Appendix I; the y, functions have been obtained by the method of variation of 


parameters, leading eventually to a numerical integration of Bessel functions of 
order 1/3, while the functions y, and y, have been obtained by inspection. 
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It is also shown that 
(12) 
Similarly ys 4 y (14) 


(15) 


Substituting equation (11) into (5), and then using equation (9), it can be shown 
that the bending moment distribution in the strut is given by 


+C,{(1 + aux) y, — px} + {(1+ y,—(uxy}. (16) 


The functions y, to y, can be described as deflection functions for a strut with a 
linearly varying axial loading P= P,,+ qx, applied along the straight strut axis, and a 
distributed lateral loading w. The deflection functions y,, vy, and y,, together with 
their first derivatives (slopes) are plotted in Figs. 3 to 8 for various values of the 
parameter z=(q°E//P,°)''®. The y, and y, functions are shown in Figs. 9 and 10, 
but for practical applications it is generally preferable to calculate these functions, 
and their first derivatives, from equations (12) to (15). In all cases the range of 
values of z has been restricted to0 <= 2<2. However, the special case of z= 00, 
corresponding to P,, =0, is discussed separately in Section 4. 

It can be seen clearly from Figs. 3, 5, 7 and 9 that as 2 tends to zero the 
deflection functions degenerate into the general solution for the deflection of a 
uniform strut subjected to constant axial loading, i.e. y,=cos ux, y,=sin ux, v,=1 
and y,= x. 


4. Deflection Functions for n=n.+qy and P,=0: Loading Type (a) 


If the end load P, =0, then «x can no longer be used as an independent variable. 
Introducing therefore Ax ={q/EI)''*x as the new independent variable, equation (8a) 
takes the form 


} w 
- +X = 7 
d(ax? d(x? *" gq’ (7) 


where A°=q/El and » is used in place of y to avoid confusion with solutions for 
P, 0, i.e. for < 00, 
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FIGURE 3 


FIGURE 4. 


FiGuRE $5. 
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FIGURE 6. 


-2t 


FIGURE 3-6. Variation of the deflection functions y and their first derivatives (slopes) y with 
px and 
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FIGURE 8. 
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Ficures 7-10. Variation of the deflection functions y and their first derivatives (slopes) y’ with 
px and gq. 
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Equation (17) can be integrated directly twice to give 


2A 


wnere C and D are arbitrary constants of integration. The complementary functions 
of this second-order differential equation are given by‘ 


[ | (19) 


where J,, and Y,,, are Bessel functions of the first and second kind respectively 
and of order 1/3. Differentiating equations (19%) and (20) with respect to (Ax) it 
follows that 


n, =dn,/d(Ax)=Ax E (Axy (21) 


and dn, /d(Ax)=Ax Y- E (Axy (22) 


As in the previous case, the method of variation of parameters can be used to 
evaluate the particular integral », when the right hand side of equation (18) is equal 
to 1 (see Appendix I); the other particular integrals can be obtained directly by 
inspection from equation (18), so that 


Hence the complete solution of equation (17) becomes 


Ax 
) = An, Bn, r Cn; Tr D rw ig (25) 


where A, B, C, and D are arbitrary constants of integration. 


Substituting equation (25) into (5), and then using equation (18), it can be shown 
that the bending moment distribution is given by 


M 


q 


Ax (An, + Bn.) +C (Axn, — 1). (26) 


The deflection functions »,, », and », and their first derivatives (slopes) are 
shown in Figs. 11 and 12. 
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Ficure 11 anp 12. Variation of the deflection functions » and their first derivatives (slopes) 
yn with Ax 


5. Deflection Functions for n=—n.: Loading Type (b) 


For n=n. equation (6) takes the form 


de El + . (27 


which becomes, using the non-dimensional variable «x, 


a? { d?y } dy w 
d(ux¥ Ud (uxp C+ d (ux) uP, 
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Equation (27a) can be integrated immediately to give 


a y 
d (uxy 


+(1+ =C,+ ux, . . 


where C, is a constant of integration. Substituting now ¢=dy/d(ux) into 
equation (28) it follows that 


d*o w 
which is of the same form as equation (9) and therefore its complete solution is 
given by 


dy wy, 


Thus in this case the deflected shape of the strut is determined by 


ut 


ax 
y=C, | y, d(ux)+C, vad (ux) +C, d(ux)+ C.4 | y,d(ux), (31) 
6 


0 


where C, to C, are constants of integration. The bending moment distribution is 
obtained from equations (5) and (31) and is given by 


M = = ( ) 

Integrals of the deflection functions y,, y, and y, have been obtained by 

numerical integration and these are shown in Figs. 13, 14, and 15. To calculate 
integrals of the y, functions, equation (12) can be used and gives 


ax mx px 
| y,d (ux)= y,d (ux)+ yd (ux).. (33) 
0 0 0 


Similarly, if the end load P,=0 and Ax is used as an independent variable, it 
can be shown that 


Ax Ax Ax 

n= A ,d(Ax)+B nd (Ax)+C (Ax)+D+w = (34) 
0 
and an = —(An,’ + Bn,’ +Cn,’), ‘ (35) 


where primes (’) denote here derivatives with respect to Ax. Integrals occurring in 
equation (34) have been obtained by numerical integration and are shown in Fig. 16. 
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Ficure 15. 
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Figures 13-15. Variation of integrals [va with ux and z. 
0 


February 1960 83 


ALC... i 
1 2 075 
-0.5 . 
-1.0 
2.0 
ors \ 
1.0 > = 
0.5 
20 
1 2 3 4 5 
azo 
3 
2 
1.0 
LS 
2.0 
Pa 3 5 6 


I S. PRZEMIENIECKI 


3 

2r 

1} 

oO A. A. | 

i 2 3 —__4 5 


Af 
Ficure 16. Variation of integrals (nd (Ax) with Ax 


6. Statically Determinate Problems 


If a strut is subjected to lateral loading or end bending moments, then it is 
possible to determine deflections, and hence bending moment distributions, for given 
end conditions. In such cases the end conditions provide a sufficient number of 
equations, at least one of which is non-homogeneous, to determine the unknown 
constants of integration in the general solution for deflections. Multi-bay struts 
present no difficulty, apart from the increase in the number of unknown constants 
of integration. 


Several common cases occurring in the analysis of aircraft structures for which 
n=n.+qy, Le. type (a), have been compiled in Appendix III. Furthermore, to 
facilitate the calculation of bending moments for the case of eccentrically loaded 
struts, which is usually the most common case in practice, the y, and y, functions 
have also been plotted in a carpet form in Figs. 17 and 18. 


7. Instability Loads 


If no external lateral loading or bending moments are applied to the strut then 
it is not possible to determine the deflections from the corresponding kinematic 
boundary conditions at the ends of the strut. The boundary conditions give then 
only a set of four homogeneous linear equations for the deflection functions y, to y, 
with a trivial solution y=0. For the deflections to have values other than zero, the 
determinant formed by the coefficients from the simultaneous equations must be 
equal to zero. When expanded, this determinant yields a transcendental equation 
which constitutes the required stability criterion for the compressive load P,, and the 
distributed loading q, and the roots of this transcendental equation give a functional 
relationship between P, and q which determines the stability boundary for the 
applied loading. 

The stability boundaries have been calculated for typical end conditions and 
are plotted in Fig. 19 as a series of curves relating the maximum compressive axial 
load P, with the maximum distributed loading gq which, when applied simul- 
taneously, produce flexural instability of the strut. This shows that struts with 
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Ficure 17. Carpet plot for the y, functions. 
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Ficure 18. Carpet plot for the y, functions. 
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Ficure 19. Stability boundaries. 


distributed loading applied along the undeformed strut axis are more stable than 
struts with loading applied at the strut axis in the deflected position; the difference 
between the two cases is a maximum when P, =0. 


8. Numerical Example 


As an example, consider a pin-ended strut subjected to bending moments 
M,=P.,e and M,=(P,+qle at x=0 and x=I/ respectively, a distributed shear 
loading g=800 Ib./in. applied along the strut axis, and axial load P,=20,000 Ib. 
Thus for this case the distributed couple n=n.+qy, where n.=qe, and e can be 
regarded as the eccentricity of applied loading. The flexural rigidity EJ = 1-673 x 
10° Ib.in.? and the strut length /=20 in. Hence 


20, 20? am Be 
s_ 800x20° _ 4. 
= joe = 3826, N=1°564, 


and a=A/u=0°715. 
Using equations (11) and (16) with w=0, the general expressions for the 
deflection and bending moment distributions become 


=(1+ uxXC,y, +C.y.)+C, [(1+ ux) y, 1] +C, (1+ 2° 4x) pr). (37) 
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These equations are to be solved subject to the boundary conditions 


y=0, M=P,e at x=0 
(38) 


y=0, M=(P,+qlDe 


Using these boundary conditions to determine the unknown constants of 
integration C, to C,, and noting from equation (12) that y,=(1—y,)/2°+y,, it 
follows finally that 


M=P, e (1+ 2 xx) [y, 1) (39) 


where the additional suffix / denotes the values of y,, y, at x=/. Variation of the 
bending moment with x//, calculated from equation (39), is shown in Fig. 20. To 
illustrate the difference between the distributed and constant loading cases the 
bending moment distributions for P= P, and P=P,,+ ql are also included in Fig. 20, 
which shows clearly that in this case the use of a constant axial load P=P,+ ql, 
i.e. P=36,000 Ib., would overestimate the maximum bending moment by a 
large factor. 


20 


6 


Q 
=z 


NON-DIMENSIONAL BENDING MOMENT 
@ 


0.4 0.6 0.8 
NON-DIMENSIONAL DISTANCE (2¢/£) 


Ficure 20. Numerical example: bending moment distribution. 
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Appendix I 
GENERAL SOLUTION OF EQUATIONS (9) AND (18) 


The governing differential equation for calculating the deflected shape of a loaded 
strut is of the form (see equation (9) ) 


where F (ux) is a function of the independent variable »x. The two complementary 
functions of this equation can be obtained by the well-known method of series integra- 
tion. It can be shown easily that this leads to 


where C, and C, are arbitrary constants of integration and y, and y, are given by 


— 3 l ( 3 4 5 9 7 
y, =cos ux+a® 4 31 + 51 (ux) — 
4 28 
+ (ux)*— (ux)* + } 
2 6 


These expressions can be used for small values of ux, but for ux > 1 the convergence 
of the series is not satisfactory. For 2=0, i.e. for constant axial loading, the two 
solutions become simply y,=cos ux and y,=sin ux. An alternative solution can be 
obtained in terms of Bessel functions of the first and second kind by changing the 
independent variable «x into a new variable ¢ given by 


Hence, provided that z is finite and not equal to zero, equation (40) with F (ux)=0 can 
be transformed into 


fy 
= =0, ; ‘ 44) 
+ 
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The solution of equation (44) is given by 
Yor=AP, (ux)+ BO, (ux), . ; (45) 


where A and B are arbitrary constants and 


2/ 
3\ a? 323 


Differentiating equations (46) and (47) with respect to «x it follows that 


and (ux)= 2( “4 ; (50) 


The series solutions y, and y, must obviously be linear combinations of the Bessel 
function solutions @, (ux) and 4, (ux), because the latter are solutions of the same 
differential equation (40). Hence 


= = (51) 
where a, b, c and d are constant coefficients. From equations (41) and (42) it follows 
immediately that 


y,(O)=1 ; y,(0)=0 | 


(52) 
y(0)=0 ; y,(O)=1 


The unknown coefficients from equation (51) can now be determined from (52) and it 
is found that 


a=—,'(0)/W 


b= /(0)/ W 


(53) 
c= ,(0)/W 
d=—®,(0)/W, 

where W= PO). 
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The expression (54) for W is the Wronskian determinant of equation (44) and it can 
be shown that it is equal to —32/x (see Appendix II). 


Bessel functions of the first and second kind of order 1/3 and —2/3, occurring in 
equations (46), (47), (49) and (50), have already been tabulated®: ® and therefore it has 
been possible to compute the y, and y, functions and their derivatives (slopes) from 
equations (51), (53) and (54). 


It remains now to determine the particular integrals of equation (40) which, 
together with the solutions y, and y,, will give four complementary functions and one 
particular integral of equation (8a). The particular integrals of equation (40) can be 
best found by means of the method of variation of parameters. This method gives the 
general solution of equation (40) in the form” 


ax 

y= | F (ux) (ux) F(ux)yd (ux) 
k k 


(55) 


where the lower limits k, and k, are arbitrary and they introduce constant multipliers 
of y, and y,, while the upper limits give the required particular integrals. It is shown 
in Appendix II that the Wronskian determinant of equations (40) or (9) is given by 
y,¥.—y,'y,=—1. Hence the calculation of the particular integrals of equation (9) 
reduces to the integration of 


Y= F (ux) y,d (ux)+y, F (ux) y,d (ux), (56) 


0 


where, for F (ux), the values 1, ux, and (ux)? are substituted in turn. For the purpose 
of computing integrals in equation (56) it has been assumed arbitrarily that k,=k,=0. 
Thus equation (56) represents the particular integral, together with some undetermined 
parts of the complementary functions y, and y,. 


Changing the independent variable x into z in equation (56) by means of equation 
(48), substituting F (ux) =1, and then using equations (46), (47), (51), (53) and (54), it 
follows that 


f f 
{ (2) J dz+J,,5 (2) | (2) de } 


2 


0 0 


with = =2/ (325), ; ; . (48a) 


where A and B are constants and y, is the particular integral of equation (40) with 
F (ux)=1. The constants A and B are determined from the condition that when 
a tends to zero, y, should tend to one, which is the correct particular integral for 
y”+y=1. Substituting y,=cos ux and y,=sin ux into equation (56) with k,=k,=0, 
integrating the resulting equation and then using equation (57), it can be shown that 
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Differentiating equation (58) with respect to ux, 


3 


z 


1 (3z\2 


z 


0 

Functions represented by equations (57) and (60) have been computed by numerical 
integration of { J, , (z)dz and {| Y, , (z) dz, since no simple expressions could be derived 
in terms of known tabulated functions. 


For F(ux)=ux and F (ux)=(ux)*, there is no need to calculate the particular 
integrals separately, since the results for F (ux)=1, i.e. the y, functions, can be directly 
used for that purpose without any additional numerical integration. The method used 
becomes apparent from Table I. 


TABLE I 
TABLE OF PARTICULAR INTEGRALS 
Equation | Differential Equation Particular 

No. Integral 

(i) y” y=1 Ys 

(ii) y” +(1+2% 4x) y=ux 
(iii) y” +(1+a% 4x) y=(uxP 
(iv) +(1+a3% y=14+ ux | 

(v) y” + (1+ 43 ux) 23 ux 


The particular integrals for equations (iv) and (v) in Table I have been obtained by 
inspection. Substituting now appropriate particular solutions into the differential 
equations (i), (ii) and (iii), and then adding equation (i) to z*® times equation (ii), 


(y,”" + ")+(1 + 22 ux) (y, + J=14+ 2% ux. 


Hence, using the particular solution for equation (iv), it follows that 


i—y, 
x3 
Comparing now the functional behaviour of this equation at the origin x=0 with the 
particular solution for «=0, i.e. y,=,x, it is found necessary to modify it by adding 
the complementary function y, so that 


+ (12) 
Zz 
, , 
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Similarly it can be shown that for F (ux)=(ux)* the particular solution y, is 
given by 


uUx—Yy, 


as 
Now from equations (9) and (8a) it follows that the functions y,, y,, y, and y, 


are the four complementary functions of equation (8a), while y, is the particular 
integral for uniformly distributed lateral loading w. Thus the complete solution of 
equation (8a) can be written as 


wy. 
y=C,y, tC.y, + C,y, tC. y, + = 2p ‘ (11) 
) 


For the case of P, =0, i.e. z= 00, the governing differential equation is 


w 
+Axn=C+DAx+ -(Ax)... (18) 
d (Ax)? 2Aq 


To evaluate the particular integrals of this equation the method of variation of para- 
meters can also be used. For the right hand side of equation (18) equal to one, the 
particular integral is given by 


Ax 
no (Ax) (Ax) 
— 12 | 
0 0 


(61) 


where the lower limits have been chosen arbitrarily. In this case the Wronskian deter- 
minant is given by »,'n,—1,'",=—3/= (see Appendix II) so that, substituting equations 


(19) and (20) into (61), it follows that 
AX 
(Ax)'/? {J, E (Axy’ | CARY | ~ 


Ax 


0 
and thence 


Ax 


Ax 
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Integrals occurring in equations (62) and (63) can be transformed into integrals 
§J,,,(6) dé and Y,,, (4) d6@ by the substitution 6=(2/3)(Ax)’/*. Since these integrals 
had to be calculated for equations (57) and (60), no additional computation is required 
for the calculation of the n, and 7,’ functions. The remaining particular integrals can 
be obtained by inspection; for the right hand side of equation (18) equal to Ax, 


(23) 
while for (Ax)? the particular integral 
(24) 


The deflection functions »,, ,, 7,: and », are the four complementary functions 


of equation (18), while 7, is the particular integral for uniformly distributed lateral 
loading w. Hence the general solution of equation (18) becomes simply 


Arn, +Bn,+Cn,+Dn,+w (25) 


~ 
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WRONSKIANS OF EQUATIONS (9), (18) AND (44) 


Let y, and y,, be two solutions of the second-order differential equation 


ly 


dx? dy 


where f (x) and g¢ (x) are functions of an independent variable x. The Wronskian W of 
the two independent solutions y, and y, is defined as 


dy, 
= ¥,. 
72 
dx dx 
This determinant can be obtained, except for a constant factor, directly from the differ- 
ential equation by means of the relationship 


z 


W=Aexp[~[ 


where A is a constant factor. Hence the Wronskians for equations (9), (18) and (44) 
are independent of x, since in all three equations f (x)=0. 


Consider now the Wronskian of equation (44), which is given by 


?,’D, 


=— J/3 az [J,,3 (2) 
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2 


and primes (’) denote derivatives with respect to ux. The last expression in equation 
(67) is obtained by means of the following identity for Bessel functions of the second 
kind : — 


{J, (z) cos vz—J_, (z)} 
sin vr 


(2)= (68) 


Since in this case the Wronskian is independent of ux and z, W can conveniently be 
calculated at z=0. Substituting therefore approximate expressions for Bessel functions 
when z is small, and then letting z tend to zero, it can be shown that 


Using now the property of the gamma function that 


sin =x’ 


equation (69) becomes 


W=—32/". . ; (71) 


Similarly the Wronskians of equations (9) and (18) can be evaluated at x=0. 


Appendix III 
BENDING MOMENTS AND DEFLECTIONS 


The formulae for the deflections and bending moments in Table II may be applied 
to any strut with constant cross section, in which the elementary formula connecting 
curvature and bending moment holds and the distributed axial loading is applied along 
the straight loading axis, i.e. loading type (a). 


The formulae for the deflections and bending moments in Table II, with the 
exception of Cases 3 and 8, have been derived directly from the basic equations (11) and 
(16) or (25) and (26), using appropriate boundary conditions. For Cases 3 and 8, where 
a concentrated lateral load Q is applied, two separate sets of basic equations with 
eight arbitrary constants must be considered, corresponding to the two regions for 
O<x</, and /,<x</. Some simplification, however, may be achieved by con- 
sidering the differential equations for the bending moment in a pin-ended strut subjected 
to a concentrated lateral load. 
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TABLE II 
SUMMARY OF FORMULAE FOR BENDING MOMENTS AND DEFLECTIONS 
(Bars over functions denote the values at x=1) 


Ends pin-ended. Compression end loads, 
distributed shear and end couples 


Mode of Support and Loading Deflection and Bending Moment 


=M,{1 + (i+ 473 


(M,—M_) { 

ul 
Ends pin-ended. Compression end loads,| {> (5 My 
distributed shear and uniform lateral, 
loading 


Real 


w 
1 


l 
+> 
ul 1 
M4 Ys , Bx ~ (ux)? 


y=Ay,+ By, P for ax <= 


Ends pin-ended. Compression end loads, . I 4 | 

distributed shear and concentrated lateral | *?¥2 
load 


Ql, gi, | 


forO< ex 


M Ql y 


0 


+ for ex <= wl, 
0 


where A, B, C and D are constants determined 
from the conditions that y=0 at each end, and | 
the conditions that y and y’ should be con- | 
tinuous at x=/.. 
Ends pin-ended. [Eccentrically applied | 
compression end loads and distributed 
shear | { | (1-5), } 
yt+e=e y,+ —=—), 
/ qf 1 2 


/ 


M=eP (1+ {», 
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TABLE II 


(continued) 


: Mode of Support and Loading 
Arbitrary end conditions. 
end loads, distributed shear, end couples, 
and uniformly distributed lateral loading 


6. Ends pin-ended. 


Ends pin-ended. Distributed shear and 


uniform lateral loading 


4 


Ends pin-ended. 
concentrated lateral load 


Deflection and Bending Moment — 
w 


M| 
py (C,y, + + + auxdy,—1) +) 
+C, x) y,— + 

w 
{(1+ y, —(ux)}, | 


where C, ....C, are arbitrary constants deter- 
mined from the boundary conditions. 


M. 
n= +0813), —Aly,)  +0°813y, — 
ql 


M, ", 
—An, +1 (1 +0°813%,)—! —0°813», i} 
ql 
M=M, { -(1 +0813, ( 
2 
x 
+M, {a 0813(*) ns} 


wil x 
0-813 —* », 08137, + ( ) 
l 


1 
M GY —n, } 
n=An, + By wae for OD 
2 gl 
Ql, 
n=Cn, + Dn, — Aln,—1) for 


Ax (An, + Bn,) for Ax 


q/A 
M Ql Ql 

Ay Cn, + Dn, — — ) 
q/A ( "2 ql ql 


for 


where A, B, C and D are constants determined 
from the conditions that »=0 at each end and 
the conditions that » and 7’ should be continu- 
ous at x=/.. 
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TABLE Il 
(continued) 


of Support and Loading 


Deflection and Bending Moment 


9. Ends pin-ended. Eccentrically | applied 
distributed shear 


af n+e=e {a+08135)™ } 


e — en 

M=eqx {a +0°813n,) -0-8131, } 
x 


10. Arbitrary end conditions. Distributed if 
shear, end couples, and uniformly distri- 
buted lateral loading 


w 
1= An, + Bn, +C1,+D+ —~ (Ax) 
q 


= Ax (An, + Bn,)+C (Axn, — 1), 


where A, B, C and D are arbitrary constants 
determined from the boundary conditions. 


For Case 3 these equations are given by 


l, 
d’y 


Writing these equations in non-dimensional form as 


y”+(1+ 2% ux) y= — ul, (76) 
and y’+(1+ ux) y= — ul, Spx < pl (77) 
| 


it becomes clear that the general solutions for the deflections are given by (see 
equation (9)) 


l 
y=Ay,+By,— AE Vas ul, : (78) 
is ulP 
l 
and y=Cy,+Dy,— al al, Spx S pl. (79) 
A 


Thus only four, instead of eight, arbitrary constants appear in the equations for deflec- 
tions. These constants are determined from the conditions that y=0 at each end and 
the conditions that y and dy/dx should be continuous at x=/,. The bending moment 
equation then follows simply from equation (5). 
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Some Properties of the Admittances 
of Dynamical Chains 


W. J. DUNCAN 


(University of Glasgow) 


Summary : Simple relations among the coefficients appearing in the admittances 
of dynamical chains are obtained by considering the initial motions which 
follow the application of certain impulses. 


1. Introduction 


This paper is a sequel to Ref. | and relations among the coefficients appearing 
in the impulsive or indicial admittances of the simple dynamical chains considered 
there are obtained from an examination of the initial motions which follow the 
application of certain single impulses or systems of impulses. The reader may be 
reminded that the dynamical chains considered consist of a set of masses (or 
generalised masses) linearly arranged and with “elastic” connections between 
adjacent elements. The masses are m,, m,...... m,, and the corresponding dis- 
placements x,, x, ....... x, respectively. The stiffness of the connection between m, 
and m,., is denoted by c,, ,,, while that of the connection between m, and “earth” is 
C,,; the stiffness of the connection between m, and earth is c,,.. When both c,, and 
Cn) are Zero the system is called unanchored: otherwise it is anchored. For anchored 
systems with n degrees of freedom the typical impulsive admittance is 


n 
p=1 


where , .... ©, are the roots of the characteristic equation. An unanchored 
system with nm degrees of freedom has the typical impulsive admittance 


n-1 
2,,(t)=A,t+ SiN o,f, 
p=1 


where A, is the reciprocal of the sum of the masses of the system, while Pinin 
are the roots of the reduced characteristic equation. The equations for determining 
the coefficients A,,, were given in the earlier paper, together with explicit expressions 
for the impulsive admittances for two- and three-mass anchored and unanchored 
systems. 
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2. Unit Impulse Applied to " Mass 


The system is initially at rest and all the displacements zero (elastic connections 
unstrained). Unit impulse is then applied to m, at the instant 1=0, so it acquires 
the initial velocity |/m,. Consequently, if we pay attention at present only to the 
lowest power of t appearing in the displacement, we shall have 


t 
(3) 
Now consider the adjacent mass ,,,. Its neighbour 7, exerts no force upon it 


initially but, on account of (3), the mass m, exerts a force upon it given by 


Hence the initial motion of m,,, is governed by the equation 


(5) 
while x,,, =x .4,=0 when t=0. Thus when is very small 
When the system is anchored the complete expression for this displacement is 


where 7 is the fumber of masses or degrees of freedom. Now expand the sines in 
powers of ¢ and compare with (6). We get 


~ 


p=1 


Cr or 
p=) m1, 


For an unanchored system with mn masses we obtain similarly 


p=! 
n Cre 
p=1 
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DYNAMICAL CHAINS 


Next consider the motion of m,,,.. The mass m,., initially exerts no force on 
m,,,, but m,,, exerts the force 


3!m,m,4, 


Cr+ 


on account of (6). Hence the expression for x,., with only the lowest power of tf 
retained is 


5!m.m,4, 


On comparison with the complete expression for x,.., namely 2,.., (0), we obtain 


the following results. 


For an Anchored System 


Cr 


M,M,. 


& 


p=1 
For an Unanchored System 


(17) 


n 


Oy Ay . (18) 
M,M, 


Evidently the argument can be extended to the (r+ k)"" mass (where r+k } n) and 
zero summations similar to (13), (14), (16) and (17) will be obtained. The first non- 
vanishing summation will be 


On account of the symmetry of the admittances we have 


A rep = Ary. 
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n 
0A. 
pat 
n—1 


W. J. DUNCAN 


By use of this relation summations containing A,,, with s less than r can be obtained 
from those already given. 


It only remains to consider the initial motion of the mass itself. Its initial 
velocity is 1/m,, so we have for an anchored system 


| 
while for an unanchored system 
5 


When ¢ is very small the force on this mass is 
(c,. + Cror4 ar 


It follows that 
m, = Wp A rrp = (c, rrp . . . (23) 


All the foregoing equations may be checked and exemplified by the results given in 
the earlier paper for two- and three-mass systems. 


For an anchored system with n masses the greatest value of k is(m—1). Then 
there are n linear equations (terminating with (19)) for the n coefficients involved in 
the admittance and these may therefore be found. Thus the admittance 


Qin (O=2,, 


is obtained completely. For an unanchored system there are only (n—1) unknown 
coefficients in each admittance, since A, is known. Hence the coefficients can be 
found from equations similar to the set (16)-(18) when k=n—2 and when k=n-I. 


3. Impulses Proportional to the Respective Masses Applied 
Simultaneously 


Impulses are applied simultaneously to all the masses in the same sense and the 
measure of the impulse applied to any mass is equal to the measure of this mass. 
Then all the masses initially have unit velocity. Moreover, when the system is 
unanchored the masses will continue to move uniformly with unit velocity and 
therefore without relative displacement. Hence, for an unanchored system, 


a 
1 
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DYNAMICAL CHAINS 


For an anchored system the condition that the initial velocity of m, is unity yields 


n n 


=> mweA,,=1. . ‘ 


p=l 


In this motion the extreme masses m, and m, are subject to accelerations initially 
proportional to ft, while the accelerations of the other masses are initially 
proportional to higher powers of t. Relations among the coefficients can be obtained 
by detailed arguments conducted along the lines of those given in Section 2, but 
these do not appear to be of much interest and will not be quoted. 


REFERENCE 


1. Duncan, W. J. Impulsive and Indicial Admittances of Simple Dynamical Chains. Aero- 
nautical Quarterly, Vol. 1X, p. 1, February 1958. 


Correspondence 


Impulsive and Indicial Admittances of Simple 
Dynamical Chains 


At the end of the paper by PRoFESSOR DuNCcAN, published in the February 1958 
issue of The Aeronautical Quarterly (Vol. LX, p. 1), the author writes that it has not so 
far been proved that all the coefficients A,,, occurring in the direct admittances z,, (t) 
are positive. A correspondent, MR. RoBerT J. FITZGERALD, of Massachusetts Institute 
of Technology, has now supplied a proof of this, as follows. 


If we express the forced-vibration mode shape [A] of a system as a sum of the 


2 [M] 02-0? 


is the set of applied force amplitudes at frequency w, 


and [M] is the inertia matrix. The superscript “T” denotes that a matrix is to be 
transposed. 
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natural mode shapes [A,] (i=1, 2,..., we get 
a, 
a, 
where [A]= [AJ=| . 
a, ay 
i, 


R. FITZGERALD 


In the case of the direct admittance, z,,=a,/f,, we let f,=f,= ... =f,=0, and 
obtain 


in which C, > 0 (i=1,2,...., n), since the denominators are positive due to the positive 
definite nature of [M]. 


If we replace w* by -s? in z,,(w), we get the Laplace transform of the impulse 
response 


as can be shown by taking Laplace transforms in the original matrix dynamic equation. 
Taking the inverse transform, the impulse response becomes 
bed 
2,,(1)= > — SiN 


in which all the coefficients are positive. 
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BRISTOL | 
Aircraft 


EXPANDING COMMITMENTS LEAD TO AN 
EXPANDING GUIDED WEAPON TEAM 
there are opportunities for— 


ADVANCED WORK ON 
HYPERSONIC PROJECTS 


AERODYNAMICS— 


Acronautical Engineers or 
Mathematicians 


Al. Advanced problems of 
flow and their application to configuration 
design. 


hypersonic 


A2. Unsteady flow werodynamics. 

A3. Stability and response. 

A4. Supersonic project aerodynamics. 

A5. Hypersonic flutter. 
STRUCTURES— 


Engineers or Mathematicians 
Sl. Project design of advanced missiles. 


82. Heat transfer problems. 


INSTRUMENTS— 
Physicists 
Pl. Optical devices, including mosaic 
systems, colour densitometry and flight 
camera installations. 


P2. Gyroscopes, pressure measuring de- 
vices and clectro-mechanical instruments. 


FLIGHT SYSTEMS— 
Fl. Mathematicians or Engineers for 
theoretical studies and development on 
missile servo-control systems, 


F2. Mechanical Engineers for theorct- 
ical studies and development on hydraulic 
and fuel systems. 
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F3. Physicists or Electronic Engineers 
for radome development. 


F4. Physicists or Electronic Engineers 
for aerial development. 

ELECTRONIC TECHNIQUES— 
Physicists 
El. Development of techniques associated 
with advanced electronic devices. 


WEAPON SYSTEMS— 
Mechanical and Electrical Engineers 


Wi. Liaison with manufacturers of the 
Bloodhound Weapon System in connection 
with overseas sales. 
These extremely interesting opportunities are 
offered by a member of Europe’s largest 
guided weapon consortium, which has 
recently received large orders for Bloodhound. 
Applicants should have a degree or similar 
qualifications and 2 or 3 years’ experience in 
related fields. Salaries will match the high 
standards required and assistance will be 
given with housing and removal expenses. 
Please write to Mr J. Raimes, giving full 
details of age, qualifications and experience and 
quote the reference numbers of the posts in which 
you are interested. Further information on these 
posts will be supplied on request. 


J. Raimes, Personne! Manager, 
Bristol Aircraft Ltd, Room AEO/60/AQ, 
Filton, Bristol. 
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V-bombers fly over 10 miles a minute, 


over 10 miles up, 
for over 10 hours.. 


... AND BRISTOL SIDDELEY 


SUPI 


One of the largest manufacturers of motive 
power units in the world, Bristol Siddeley 
Engines Limited produce two outstanding 
high-thrust turbojet engines—the Olympus 
and the Sapphire. Between them these two 
engines power the major part of the RAF’s 
strategic V-bomber force. The Olympus, 
which delivers 17,000 Ib thrust dry (other 
versions are rated at 33,000 lb), powers the 
Avro Vulcan Mk 1 and 2. The Sapphire 
powers the Handley Page Victor Mk 1. 
These engines give the V-bombers sonic 
capability—long range—great altitude— 
superior performance to any other aircraft 
of their type in the world. 
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LY THE POWER 


Bristol Siddeley products currently in produc- 
tion include : 


AERO 
TURBOJETS: Olympus, Sapphire, Orpheus, 
Viper. 
TURBOPROPs: Proteus, Double Mamba. 
ROCKETS: Gamma. 
RAMJETS: Thor. 


INDUSTRIAL & MARINE 
DIESELS: Maybach. 
GAS TURBINES: Proteus. 
cars: Armstrong Siddeley Star Sapphire, 
Bristol 406. 
OTHER PRODUCTS: Gas bearing compressors, 
Beaver ball screws and splines. 


BRISTOL SIDDELEY ENGINES LIMITED 
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‘50 years of experionce 
are, behind the Blackburn NA 


if you are an A.F.R.Ae.8. or possess an equivalent 
qualification, please write to the Technical Staff Manager, 


when we shall be pleased to give you details of 
the opportunities available to experienced men within our 


technical organisation. 


Blackburn Aircraft Limited, Brough, Yorkshire 


Al28/d 
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aircraft at Fim Yorkshire. 
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